import os
import sys
import math
from datetime import date,datetime,time
from xlrd import open_workbook,xldate_as_tuple
import numpy
import scipy
from scipy.optimize import curve_fit, minimize
from scipy.stats import linregress
import cPickle as pickle
import matplotlib
matplotlib.use('agg')
from matplotlib import pyplot as PLT
from copy import deepcopy

"""
This module contains classes and functions designed for the automated analysis of CalTorque data.

Classes:
    CalTorque(excelFilename): create single run data object from excel file
    RunHistory(): save run history information to pickle
    SourceTube(): generate source tube object

Functions:
    OpenCalTorqueDB(path): open database
    SaveCalTorqueDB(data, path): save to database
    GetAllRunsFromFolder(path): save all runs in folder to database
    GetNewRunsFromFolder(data,path): save runs in folder that are not in data to database
    PlotCalTorque(run, tube, outdir): plot torque vs time and torque vs position
    PlotVelocity(run, outdir): plot velocity vs time
    PlotAvgTorque(data, average, dir, bins, max_d, outdir): plot torque vs position data average of all runs
    PlotSources(data, src, bins, max_d, outdir): plot average torque vs position for each source
    PlotMonthlyAvgs(data, months, bins, max_d, outdir): plot average torque vs position for each month
    PlotTelescope(data, telescope, bins, max_d, outdir): plot average torque vs position before and after Feb22, 2012
    PlotMaxtvsDate(data, outdir): plot maximum torque vs date of deployment for each run
    PlotMaxtvsVel(data, bins, max_d, outdir): plot maximum torque vs velocity for each run
    PlotOffsetvsP0(data, outdir): plot offset vs initial position for each run
    PlotParameters(data, outdir): for each Poly2Gauss parameter, plot parameter vs date for each run
    PlotFit(run, outdir): plot Poly2Gauss fit to position vs time data
    
"""

# EXCEL uses very weird time units, this is a second in those units
excel_timeunit = 1.1574074074079903e-05 # sec/unit

# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
#                                             CALTORQUE CLASS
# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
class CalTorque:
    
    """Class represents a single CalTorque run data object.
    
    Public methods:
        GetBinData(dl): get bin position, torque data
    
    Mutators:
        ChangeSource(): modify source
        ChangeStation(): modify station
    
    Accessors:
        Datetime(): get datetime
        Position(): get position data
        NewPosition(): get position data minus offset found by Poly2Gauss
        Time(): get time data
        Torque(): get torque data
        Velocity(): get velocity data
        IntegratedVelocity(): get integrated velocity data
        Filename(): get excel filename
        NewParms(): getPoly2Gauss fit parameters
        Offset(): get offset computed by Poly2Gauss
        Station(): get station deployed to
        Source(): get source used in run
        Direction(): get direction of run (in/out)
        FinalPosition(): get final position
        MaxTorque(): get max torque during run
    
    Instance Variables:
        self.datetime (object): Datetime object representing date and time of run
        self.dt (float, array): time difference between start of run to present data read
        self.initial_angle (int): initial angle value of run
        self.position (float, array): position data
        self.torque (float, array): torque data
        self.rpm (float, array): angular speed of deployment data
        self.filename (str): name of excel file from which data is taken
        self.direction (str): direction of deployment (in or out)
        self.station (str): station deployed to during run
        self.source (str): calibration source used during run
        self.bin_x (float, array): binned position data
        self.bin_y (float, array): binned torque data
        self.newParms (float, list): parameters of Poly2Gauss fit of position vs time data
        self.offset (float): position difference between data and model, found by Poly2Gauss method
        self.velocity (float, list): velocity of deployment data computed from rpm data. Negative indicates retraction
        
    """
    
    def __init__(self, filename):
        """Get data from excel file.
            
        Parameters:
        filename (str) -- path to and name of excel file containing CalTorque data
        """
        # Import file from excel
        self._ImportCalTorqueData(filename)
        self.filename = filename.split('/')[-1]
        
        #Determine direction
        if self.position[0] < self.position[-1]:
            self.direction = "in"
        else:
            self.direction = "out"
        
        #Get source/station for each run
        self.station, self.source = 'unknown', 'unknown'
        self._GetSourceStation()
        
        # Produced bin_x and bin_y with default binning
        self.bin_x, self.bin_y = self.GetBinData()
        
        # Get Gaussian parameters and offset
        self.newParms = None
        self.offset = 0
        if self.direction == "in" and self.position[0] < 100 and self.source != 'weakRa-226':
            try:
                self._FitPolyTwoGauss()
            except:
                pass
        
    def _ImportCalTorqueData(self,filename):
        """Import and convert data from excel file.
            
        Parameters:
            filename (str) -- name of excel file containing CalTorque data
        """
        book = open_workbook(filename, logfile=open(os.devnull, 'w'))
        # We're only interisted in the first sheet
        sheet = book.sheet_by_index(0)
        data = {}
        
        # Get column headers from row0. Expect row0 to be xlrd.XL_CELL_TEXT type
        col0 = sheet.row_values(0)
        for i,v in enumerate(col0):
            data[v] = sheet.col_values(i,1)
        
        # Get the datetime for the first data point
        xlsdate = xldate_as_tuple(data['Date'][0],book.datemode)
        xlstime = xldate_as_tuple(data['Time'][0],book.datemode)
        data['datetime'] = datetime(xlsdate[0],xlsdate[1],xlsdate[2],xlstime[3],xlstime[4],xlstime[5]) #2013-05-29 10:37:45
        
        # Compute the fractional second of elasped time
        data['dt'] = []
        for i in range(len(data['Time'])):
            data['dt'].append((data['Time'][i] - data['Time'][0]) / excel_timeunit)
        
        # Translate values to useful units and measurements
        self.datetime = data['datetime']
        self.dt = numpy.array(data['dt'],dtype='float32')
        self.initial_angle = data['Angle'][0]
        self.position = (numpy.array(data['Angle'],dtype='float32'))/360.*3.14159*4.0*2.54
        self.torque   = numpy.array(data['Tracking Value'],dtype='float32')*141.612
        self.rpm = numpy.array(data['RPM'],dtype='float32')
        SPROCKET_RADIUS = 2 #in
        self.velocity = (2 * self.rpm * scipy.pi / 60.) * (SPROCKET_RADIUS * 2.54)
        
    def GetBinData(self,dl=2):
        """Bin data to 2 cm intervals. Bin's value is max within interval.
        
        Parameters:
            dl (int) -- width of bin
        
        Returns:
            binned x and y data
        """
        
        pos = []
        for p in self.position:
            pos.append(p - self.position[0])
        xmin = min(pos)
        xmax = max(pos)
        
        bin_x = numpy.arange(xmin,xmax,dl)
        bin_y = numpy.zeros(bin_x.size)
        for i, x in enumerate(bin_x):
            tmp_data = []
            for j,p in enumerate(pos):
                if p >= x and p < x+dl:
                    tmp_data.append(self.torque[j])
            if len(tmp_data) > 0:
                bin_y[i] = numpy.max(tmp_data)
        
        return bin_x, bin_y

    def _FitPolyTwoGauss(self):
        """Fit Poly2Gauss model to data to determine the offset"""
        
        if self.station == 'S2':
            xtrn = self.bin_x[140:170]
            ytrn = self.bin_y[140:170]
        else:
            xtrn = self.bin_x[140:-10]
            ytrn = self.bin_y[140:-10]
        
        # perform poly_fit to get initial guess on baseline
        baseline = curve_fit(self._PolyModel, xtrn, ytrn)
        baseParms = baseline[0]
        
        # perform fit
        fitParms = [5, 315, 5, 5, 5, baseParms[0], baseParms[1], baseParms[2], baseParms[3]]
        fit_bounds = [(2, 20), (280, 330), (1, 10), (2, 20), (1, 10), (None, None), (None,None), (None, None), (None, None)]
        result = minimize(self._GetP2GaussResidual,fitParms, args = (xtrn, ytrn), method='TNC', bounds=fit_bounds,tol=0.0000005,options={'disp':False})
        self.newParms = result.x
        self.offset = self.newParms[1] - fitParms[1] + self.position[0]
        
    def _PolyModel(self, x, x0, a, b, c):
        """Compute second degree polynomial.
        
        Parameters:
            x (float) -- position variable
            x0 (float) -- initial position
            a (float) -- second degree coefficient
            b (float) -- first degree coefficient
            c (float) -- zero degree coefficient
        
        Returns:
            (float) Value of polynomial
        """
        return a*(x-x0)**2 + b*(x-x0) + c
        
    def _GetP2GaussResidual(self, parms, x, y):
        """Compute residual of Poly2Gauss fit with data.
        
        Parameters:
            parms (float list) -- parameters for fit
            x (float array) -- position data
            y (float array) -- torque data
        """
        model_y = self._PolyTwoGaussModel(x,parms)
        return numpy.sum((model_y - y)**2)

    def _PolyTwoGaussModel(self, x, p):
        """Compute y value for Poly2Gauss model.
        
        Parameters:
        x -- position variable
        p -- fit parameters,
             p[0], p[1], p[2] = amplitude , center_x, sigma: for peak A
             p[3], p[4] = amplitude, sigma: for peak B
             p[5] = polynomial offset
             p[6], p[7], p[8] = ax^2 + bx + c
        
        Returns:
            (float) Value of model
        
        """
        # offset = distance between peaks
        BEND_OFFSET = 52.0 #cm distance between bends gaussians are fit to
        y = numpy.zeros(len(x))
        for i, j in enumerate(x):
            y[i] = p[0]*math.exp(-0.5*((j-p[1])/p[2])**2.) + p[3]*math.exp(-0.5*((j-(p[1]+ BEND_OFFSET))/p[4])**2.) + p[6]*(j - p[5])**2 + p[7]*(j - p[5]) + p[8]
            if y[i] <= 1e-5:
                y[i] = 1e-5
        return y
    
    def _GetSourceStation(self):
        """Get source and station for a run from RunHistory pickle."""
        runHistory = OpenCalTorqueDB("RunHistory.pkl")
        for run in runHistory:
            if runHistory[run]['start_datetime'].strftime("%Y%m%d") == self.datetime.strftime("%Y%m%d"):
                self.station = runHistory[run]['pos']
                self.source = runHistory[run]['source']
                if self.direction == 'in':
                    pos = self.position[-1]
                else:
                    pos = self.position[0]
                if self._StationCheck(self.station, pos) and self._SourceCheck():
                    break
        
    def _StationCheck(self, station, finalPos):
        """Check that station is correct. For the cases where there are multiple runs in a day.
        
        Parameters:
            station (str) -- station deployed to (from RunHistory pickle)
            finalPos (float) -- final position of run
        
        Returns:
            (boolean) False if finalPos is not in range surrounding station, True otherwise.
        """
        
        if 'S2' in station and finalPos > 320 and finalPos <= 366:
            return True
        elif 'S5' in station and finalPos > 366 and finalPos <= 419:
            return True
        elif 'S8' in station and finalPos > 419 and finalPos <= 474:
            return True
        elif 'S11' in station and finalPos > 474 and finalPos <= 525:
            return True
        elif 'S14' in station and finalPos > 525 and finalPos <= 577:
            return True
        elif 'S17' in station and finalPos > 577:
            return True
        return False
    
    def _SourceCheck(self):
        """Check that source is correct. For cases when multiple runs in a day"""
        for s in ['Co', 'Cs', 'Ra', 'sTh']:
            if s in self.filename:
                if s == 'Co' and self.source == 'weakCo-60':
                    return True
                elif s == 'Cs' and self.source == 'weakCs-137':
                    return True
                elif s == 'Ra' and self.source == 'weakRa-226':
                    return True
                elif s == 'sTh' and self.source == 'strongTh-228':
                    return True
                else:
                    # if source does not match filename, it is wrong
                    return False
        # if source not in filename, assume to be correct
        return True
    
    
    def _IntegVel(self):
        """Compute velocity by taking derivative of spline fit to position vs time data.
        
        Returns:
            yprime (float array) -- array of y spline values
        """
        timeChop = []
        posChop = []
        for i, p in enumerate(self.position):
            if p == self.pos[0] or p == self.position[-1] or p >= 385: continue
            posChop.append(p)
            timeChop.append(self.time[i])
        spl = UnivariateSpline(timeChop, posChop) # s=2 works
        ys = spl(timeChop)
        d = spl.derivative()
        yprime = d(timeChop)
        return yprime
        
        
    #____________________ MUTATORS _______________________#
    # If an item is modified in the data base, the database must be re-serialised
    # Example:
    # from CalTorque import OpenCalTorqueDB, SaveCalTorqueDB
    # database = OpenCalTorqueDB("database.pkl")
    # run = database['data']['Live Graph 01-Apr-2013_in.xls']
    # run.ChangeStation("S2")
    # SaveCalTorqueDB(database)
    
    def ChangeStation(self, s):
        """ Modify station deployed to.
        
        Parameters: 
            s (str) -- new station (S2, S5, S8, S11,...)
        """
        self.station = s
        
    def ChangeSource(self, s):
        """ Modify source used in deployment.
            
        Parameters: 
            s (str) -- new source (weakTh-228, strongTh-228, ...)
        """
        self.source = s
    
    def ChangeDirection(self, d):
        """Modify direction of run
        
        Parameters:
            d (str) -- new direction ('in' or 'out')
        """
        self.direction = d
    #____________________ ACCESSORS _______________________#

    def Datetime(self):
        """Return datetime (e.g. 2013-05-29 10:37:45)."""
        return deepcopy(self.datetime)
        
    def BinData(self):
        """Return binned position and torque data."""
        return deepcopy(self.bin_x), deepcopy(self.bin_y)
    
    def Position(self):
        """Return position array in cm."""
        return deepcopy(self.position)
    
    def NewPosition(self):
        """Return position minus offset."""
        return deepcopy(self.position - self.offset)
    
    def Time(self):
        """Return time array in seconds."""
        return deepcopy(self.dt)
        
    def Torque(self):
        """Return torque readings array in in.oz."""
        return deepcopy(self.torque)
    
    def AngularVelocity(self):
        """Return angular velocity readings in rpm"""
        return deepcopy(self.rpm)
        
    def Velocity(self):
        """Return velocity of deployment in cm/s."""
        return deepcopy(self.velocity)
    
    def IntegratedVelocity(self):
        """Return integrated velocity."""
        return self._IntegVel()
    
    def Filename(self):
        """Return filename string."""
        return self.filename
    
    def NewParms(self):
        """Return parameters of gaussian fit."""
        return deepcopy(self.newParms)
    
    def Offset(self):
        """Return offset between data and gaussian in cm."""
        return self.offset
    
    def Station(self):
        """Return station."""
        return self.station
    
    def Source(self):
        """Return source."""
        return self.source
        
    def Direction(self):
        """Return direction."""
        return self.direction
    
    def FinalPosition(self):
        """Return final position from parked."""
        return self.position[-1] - self.offset
    
    def MaxTorque(self):
        """Return max torque reading."""
        return self.torque.max()
    
    def __str__(self):
        return "CalTorque[Date: " + self.datetime.strftime('%Y%m%d') + ", Direction: " + self.direction + ", Source: " + self.source + ", Station: " + self.station + "]"

# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
#                                                 RUN HISTORY
# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
class RunHistory:
    """ Get run history from xml files from local pickle (or online if pkl does not exist), 
    example entry
        <run>5872</run>
        <quality>NOT-SET</quality>
        <runType>Data-Source calibration</runType>
        <startTime>2013-11-03T15:53:45.000-0800</startTime>
        <endTime>2013-11-03T17:59:09.000-0800</endTime>
        <exposure>7521278143000</exposure>
        <comment>wTh-228 at S5, +4417deg</comment>
        <sourceType>Th-228:weak</sourceType>
        <sourcePosition>S5: P4_px ( 25.4, 0.0, 0.0)</sourcePosition>
    """
    def __init__(self):
        """Load or create RunHistory pickle"""
        #Assume there is a pickle file w/ the RunHistory. If not make one.
        try:
            self.run = pickle.load(open('RunHistory.pkl','r'))
        except:
            self.UpdateRunHistory()
        
    def UpdateRunHistory(self):
        """Get updated RunHistory from the RunHistory database online and save it locally."""
        import urllib
        import json
        openUrl = urllib.urlopen('http://exo-data.slac.stanford.edu/ExoDatacat/rest/runs?runMin=2000&runMax=')
        dataLines = json.loads(openUrl.readlines()[0])
        self.run = {}
        for i in dataLines:
            if i['runType'] != 'Data-Source calibration':
                continue
            self.run[i['run']] = i
            self.run[i['run']]['start_datetime'] = self._ParseDateStr(i['startTime'])
            self.run[i['run']]['end_datetime'] = self._ParseDateStr(i['endTime'])
            self.run[i['run']]['pos'] = self._ParsePosition(self.run[i['run']])#i['sourcePosition'].split(':')[0]
            self.run[i['run']]['source'] = i['sourceType'].split(':')[1] + i['sourceType'].split(':')[0]
        pickle.dump(self.run,open('RunHistory.pkl','w'))

    def _ParsePosition(self,run):
        """Get source position from RunHistory database
        
        Parameters:
            run (dict): RunHistory entry
        
        Returns:
            station
        """
        positions = ['S2','S5','S8','S11','S14','S17', 'Other']
        for pos in positions:
            if run['comment'] == None: continue
            # if station in comment and comment relevant to source position, get station from comment
            if pos + '+' in run['comment'] and "S5+1" not in run['comment']: # skip S5+1", as it is new S5
                return self._CheckComment(pos + '+', run['comment'])
            elif pos + '-' in run['comment']:
                return self._CheckComment(pos + '-', run['comment'])
        # if position not in comments, get it from sourcePosition entry
        return run['sourcePosition'].split(':')[0]
    
    def _CheckComment(self, pos, comments):
        """Check comment from RunHistory database for station.
        
        Parameters:
            pos (str): station in 'sourcePosition' section of database
            comments (str): comments from 'comment section of database
            
        Returns:
            station in comments, if found.
        """
        # find greatest index of station substring in comment
        start = comments.find(pos)
        # look for units - to know where to truncate comment
        for u in ['cm', '"', 'inch', "''"]:
            if u not in comments: continue
            end = comments.find(u, start)
            return comments[start: end + len(u)]
        # if no units are found, truncate comment at next whitespace
        end = comments.find(' ', start)
        return comments[start: end]
    
    def _ParseDateStr(self,string):
        """Get run datetime object.
        
        Parameters:
            string (str): date and time representation. String type: '%Y-%m-%dT% 2011-07-19T19:49:32.000-0700
        
        Returns:
            date and time as datetime object
        """
        # Cant use %z till python3.2 :: See http://bugs.python.org/issue6641
        string_datetime = string[0:-5]
        #utc_offset_sign = string[-5:-4]
        #utc_offset_hour = string[-4:-2]
        #utc_offset_min  = string[-2:]
        return datetime.strptime(string_datetime,'%Y-%m-%dT%H:%M:%S.%f')

# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
#                                            SOURCE TUBE DEFINITION
# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
class SourceTube:
    
    """Represents guide tube. Has angle and position data.
    
    Public functions:
        GetSourceTubeDC(): return deep copy of angle vs position data representing tube
        
    Instance variables:
        self.tube (dict): contains angle and position data
    """
    
    def __init__(self):
        """Generate position vs angle data representing the guide tube."""
        self.tube = {}
        self.tube['x'] = [0]
        self.tube['y'] = [720]  #2 turns * 360 = 720
        self._AddStraight(47.8)
        self._AddBendTube(29.106,90)
        self._AddStraight(136.4-93.5)
        self._AddBendTube(5.2388,90)
        self._AddStraight(196.5-144.6)
        self._AddBendTube(5.238,235)
        self._AddStraight(228.3-218.0)
        self._AddBendTube(5.2387,145)
        self._AddStraight(274.6-241.5)
        self._AddBendTube(7.618,75)
        self._AddStraight(309.0-284.6)
        self._AddBendTube(7.618,75)
        self._AddStraight(358.8-318.9)
        self._AddBendTube(7.619,90)
        self._AddStraight(414.5-370.8)
        self._AddBendTube(7.619,90)
        self._AddStraight(445.7-426.5)
        self._AddBendTube(5.08,108)
        self._AddStraight(467.5-455.3)
        self._AddBendTube(7.619377,90)
        self._AddStraight(523.8-479.5)
        self._AddBendTube(7.6193,90)
        self._AddStraight(571.6-535.8)
        self._AddBendTube(7.619,90)
        self._AddStraight(627.9-583.6)
        self._SubTailLoops(2)
        self.x = numpy.array(self.tube['x'])
        self.y = numpy.array(self.tube['y'])

    def _AddStraight(self,dx,nx=100):
        """Generate angle vs position data for straight section of guide tube.
        
        Parameters:
            dx (float): length of straight section of tube
            nx (int): number of desired data points in straight section of guide tube
        """
        x0 = self.tube['x'][-1]
        y0 = self.tube['y'][-1]
        for i in range(nx):
            self.tube['x'].append(x0 + dx*float(i)/float(nx-1))
            self.tube['y'].append(y0)

    def _AddBendTube(self,r,theta,nx=100):
        """Generate angle vs position data for curved section of guide tube.
        
        Parameters:
            r (float): radius of curved section of tube
            theta (float): angle of arc in degrees
            nx (int): number of desired data points in straight section of guide tube
        """
        x0 = self.tube['x'][-1]
        y0 = self.tube['y'][-1]
        for i in range(nx):
            self.tube['x'].append(x0 + 2.0*3.14159*r*theta/360.*float(i)/float(nx-1))
            self.tube['y'].append(y0 + theta*float(i)/float(nx-1))

    def _SubTailLoops(self,loops):
        """Subtract tail integrated bend as cable is deployed.
        
        Parameters:
            loops (int): number of loops
        """
        bend = 360.0*loops*self.tube['x'][-1]/750.0
        for i,x in enumerate(self.tube['x']):
            self.tube['y'][i] = self.tube['y'][i] - bend*x/self.tube['x'][-1]
    
    def GetSourceTubeDC(self):
        """Return deepcopy angle vs position data of Source tube object."""
        return deepcopy(self.tube)
        

# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
#                                            GENERAL USE FUNCTIONS
# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
def OpenCalTorqueDB(path):
    """Open pickle of CalTorque runs.
    
    Parameters:
        path (str): path to pickle (including filename and extension)
    
    Returns:
        pickle open for read
    """
    import pickle
    if path.split('.')[-1] == 'pklz':
        import gzip
        return pickle.load(gzip.open(path,'rb'))
    elif path.split('.')[-1] == 'pkl':
        return pickle.load(open(path,'r'))
    else:
        print "OpenCalTorqueDB: Unknown pickle extension %s" %path

def SaveCalTorqueDB(data,path):
    """Save iterable of CalTorque runs to a pickle.
    
    Parameters:
        data (CalTorque, dict): dictionary of CalTorque objects to save to database
        path (str): path to directory containing new data files
    """
    import pickle
    try:
        if path.split('.')[-1] == 'pklz':
            import gzip
            pickle.dump(data,gzip.open(path,'wb'))
        elif path.split('.')[-1] == 'pkl':
            pickle.dump(data,open(path,'w'))
        else:
            print "You did not specify a valid file extension: [pkl,pklz]"
    except:
        print "You did not specify a valid file."


def GetAllRunsFromFolder(path):
    """Add all files from path to data
    
    Parameters:
        path (str): path to directory containing data files
    
    Returns:
        data (CalTorque, dict): dictionary of CalTorque objects
    """
    #Make a consistant path 
    if path[-1] == '/': path = path[:-1]
    list_dir = os.listdir(path)
    
    data = []
    for file in list_dir:
        if file.split('.')[-1] not in ['xls','xlsx']:
            continue
        print "Importing: %s"%file
        data.append(CalTorque(path +'/'+ file))
    return data

def GetNewRunsFromFolder(data,path):
    """Process and add to data files not found in database.
    
    Parameters:
        data (CalTorque, dict): dictionary of CalTorque objects from database
        path (str): path to directory containing new data files
    
    Returns:
        data (CalTorque, dict): data dictionary updated with the new files
        newFiles (str, list): list of new files 
    """
    if path[-1] == '/': path = path[:-1]
    list_dir = os.listdir(path)
    if type(data) != dict: print "Data must be type dict"
    files_in_data = []
    for data_file in data:
        files_in_data.append(data_file.split('/')[-1])
    newFiles = []
    for file in list_dir:
        # Check extension
        if file.split('.')[-1] not in ['xls','xlsx']:
            continue
        # Check if we already have this file in data    
        if file not in files_in_data:
            print "Importing: %s"%file
            data[file] = (CalTorque(path +'/'+ file))
            newFiles.append(file)
    return data, newFiles
    
def MakeChangeCalDB(runKey, changeType, change, addToList = True, database = "CalDatabase.pkl"):
    """Method to facilitate manual changes in the CalDatabase
    
    Parameters:
        runKey (str): database dictionary reference to run (excel filename)
        changeType (int): type of data to change
            
            Key     Change
             0       change source
             1       change station
             2       change direction
        
        change (see type below): value of change to make
            
            Key     Type    Form
             0       str     '<strength><source>', e.g. 'weakTh-228'
             1       str     'S<#>', e.g. 'S5'
             2       str     'in' or 'out'
        
        addToList (bool): If True, adds this change to a list of manual changes to make to any new
                          database, default True
        database (str): path to and name of pickled database dictionary, default "CalDatabase.pkl"
    
    Note: this method permanently modifies the database. A log is kept of the changes (change_log.txt)
    """
    # open database
    try:
        db = OpenCalTorqueDB(database)
    except Exception, e:
        print str(e)
        print "ERROR: Database <%s> could not be openned."%database
        sys.exit(0)
    
    # get run to modify
    try:
        run = db['data'][runKey]
    except Exception, e:
        print str(e)
        print "ERROR: runKey did not match any in the database. Please verify that the key is correct and try again."
        sys.exit(0)
    
    previous = None
    
    # make requested change
    if changeType == 0:
        previous = run.Source()
        run.ChangeSource(change)
    elif changeType == 1:
        previous = run.Station()
        run.ChangeStation(change)
    elif changeType == 2:
        previous = run.Direction()
        run.ChangeDirection(change) # TODO add more changeable things?
    else:
        print "ERROR: changeType must be an integer in range 0 to 1. See method string for more information."
        sys.exit(0)
    
    # save database
    SaveCalTorqueDB(db, "CalDatabase.pkl")
    
    # write change made to log
    log = open("change_log.txt", 'a')
    log.write("Change on " + datetime.now().strftime('%Y%m%d-%H:%M:%S') + "\n")
    log.write("\tRun: %s\n"%runKey)
    log.write("\tChanged %s to %s in %s\n\n"%(previous, change, database))
    log.close()
    
    # add change to script to be used on new DB
    if addToList:
        fix = open("CalManualChanges.py", 'a')
        fix.write("MakeChangeCalDB('%s', %d, '%s', False, '%s')\n\n"%(runKey, changeType, change, database))
        fix.close()
    
# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--
#                                            PLOTTING FUNCTIONS
# --ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo----ooo000OOO000ooo--

def PlotCalTorque(run, tube, outdir):
    """Plot torque vs time and torque vs posistion for an individual run.
    
    Parameters:
        run (CalTorque): single CalTorque object to plot
        tube (SourceTube): SourceTube object
    """
    # get data
    position = run.NewPosition()
    torque = run.Torque()
    direction = run.Direction()
    station = run.Station()
    source = run.Source()
    time = run.Time()
    date = run.Datetime().strftime('%Y%m%d')
    # get part of filename for plot filename
    try:
        filedirection = run.Filename().split(' ')[-1].split('.')[0].split('_')[1]
    except:
        filedirection = run.Filename().split(' ')[-1].split('.')[0]
    
    # ensure folder for plots exists
    if not os.path.exists(outdir):
        os.makedirs(outdir)
    
    # plot
    fig = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax1 = fig.add_axes(box)
    ax1.set_ylabel('Position (cm)',color='b', fontsize = 18)
    ax2 = ax1.twinx()
    ax2.set_ylabel('Torque (in.oz)', color='g', fontsize = 18)
    ax1.set_xlabel('Time (s)', fontsize = 18)
    maxx = max(time)
    minx = max(time)*0.07
    ax1.hlines(338.0,minx,maxx,color='0.82',linewidth=4)
    ax1.annotate('S2',xy=(minx*.7,338.0),color='0.75')
    if max(position)> 350:
       ax1.hlines(388.0,minx,maxx,color='0.82',linewidth=4)
       ax1.annotate('S5',xy=(minx*.7,388.0),color='0.75')
    if max(position)> 400:
       ax1.hlines(442.0,minx,maxx,color='0.82',linewidth=4)
       ax1.annotate('S8',xy=(minx*.7,442.0),color='0.75')
    if max(position)> 450:
       ax1.hlines(497.0,minx,maxx,color='0.82',linewidth=4)
       ax1.annotate('S11',xy=(minx*.7,497.0),color='0.75')
    if max(position)> 530:
       ax1.hlines(548.0,minx,maxx,color='0.82',linewidth=4)
       ax1.annotate('S14',xy=(minx*.7,548.0),color='0.75')
    PLT.title(date + ' Source ' + source + ' ' + direction + ' ' + station + ' Torque vs Time', fontsize = 20, y = 1.02)
    ax1.grid()
    ax1.tick_params(labelsize=18)
    ax2.tick_params(labelsize=18)
    ax1.plot(time, position, '-', color='blue',linewidth=2)
    ax2.plot(time, torque, '-', color='green',linewidth=2)
    PLT.savefig(outdir + date + run.Filename().split('-20')[-1].split('.')[0][2:] + '_time.png')
    PLT.close(fig)

    # -- SECOND PLOT --
    fig2 = PLT.figure(figsize=(15,8),dpi=150)
    ax21 = fig2.add_axes(box)
    ax21.set_ylabel('Bend Angle',color='b', fontsize = 18)
    ax22 = ax21.twinx()
    ax22.set_ylabel('Torque (in.oz)', color='g', fontsize = 18)
    ax21.set_xlabel('Distance (cm)', fontsize = 18)
    miny = math.floor(min(torque) / 5) * 5
    ymax = max(torque)*0.5
    maxx = max(position)
   
    ## Add station locations
    if maxx >= 328.9:
         ax22.axvline(338.9, color = 'r', linewidth = 1, linestyle = '-')
         ax22.annotate('S2', xy = (340.9, miny * 0.8), color = 'r')
    if maxx >= 382.6:
         ax22.axvline(392.6, color = 'r', linewidth = 1, linestyle = '-')
         ax22.annotate('S5', xy = (394.6, miny * 0.8), color = 'r')
    if maxx >= 435.7:
         ax22.axvline(445.7, color = 'r', linewidth = 1, linestyle = '-')
         ax22.annotate('S8', xy = (447.7, miny * 0.8), color = 'r')
    if maxx >= 491.7:
         ax22.axvline(501.7, color = 'r', linewidth = 1, linestyle = '-')
         ax22.annotate('S11', xy = (503.7, miny * 0.8), color = 'r')
    if maxx >= 595.8:
         ax22.axvline(605.8, color = 'r', linewidth = 1, linestyle = '-')
         ax22.annotate('S17', xy = (607.8, miny * 0.8), color = 'r')
 
    ## Add bend locations
    if maxx >= 47.8:
         ax22.axvspan(47.8, 93.5, color = 'g', alpha = 0.2)
         ax22.annotate('Bellows', xy = (48, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 136.4:
         ax22.axvspan(136.4, 144.6, color = 'g', alpha = 0.2)
         ax22.annotate('First outer bend', xy = (136, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 196.5:
         ax22.axvspan(196.5, 218.0, color = 'g', alpha = 0.2)
         ax22.annotate('The "S"', xy = (222, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 228.3:
         ax22.axvspan(228.3, 241.5, color = 'g', alpha = 0.2)
    if maxx >= 274.6:
         ax22.axvspan(274.6, 284.6, color = 'g', alpha = 0.2)
         ax22.annotate('Bend A5', xy = (274, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 309.0:
         ax22.axvspan(309.0, 318.9, color = 'g', alpha = 0.2)
         ax22.annotate('Bend A6', xy = (309, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 358.8:
         ax22.axvspan(358.8, 370.8, color = 'g', alpha = 0.2)
         ax22.annotate('Bend A7', xy = (359, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 414.5:
         ax22.axvspan(414.5, 426.5, color = 'g', alpha = 0.2)
         ax22.annotate('Bend A8', xy = (414, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 445.3:
         ax22.axvspan(445.3, 455.3, color = 'g', alpha = 0.2)
         ax22.annotate('Bend A9', xy = (445, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 467.5:
         ax22.axvspan(467.5, 479.5, color = 'g', alpha = 0.2)
         ax22.annotate('Bend A10', xy = (467, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 523.8:
        ax22.axvspan(523.8, 535.8, color = 'g', alpha = 0.2)
        ax22.annotate('Bend A11', xy = (523, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 571.6:
        ax22.axvspan(571.6, 583.6, color = 'g', alpha = 0.2)
        ax22.annotate('Bend A12', xy = (571, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 627.9:
        ax22.axvspan(627.9, 639.9, color = 'g', alpha = 0.2)
        ax22.annotate('Bend A13', xy = (627, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 652.7:
        ax22.axvspan(652.7, 660.7, color = 'g', alpha = 0.2)
        ax22.annotate('Bend A14', xy = (652, miny), color = 'g', rotation = 'vertical', va = 'bottom')
    if maxx >= 678.5:
        ax22.axvspan(678.5, 681.1, color = 'g', alpha = 0.2)
        ax22.annotate('Bend A15', xy = (678, miny), color = 'g', rotation = 'vertical', va = 'bottom')
        
    PLT.title(date + ' Source ' + source + ' ' + direction + ' ' + station + ' Torque vs Position', fontsize = 20, y = 1.02)
    ax21.grid()
    ax21.tick_params(labelsize=18)
    ax22.tick_params(labelsize=18)
    ax22.plot(position, torque, '-', color='g',linewidth=2)
    ax21.plot(tube['x'],tube['y'], '-', color='b',linewidth=2)
    fig2.savefig(outdir + date + run.Filename().split('-20')[-1].split('.')[0][2:] + '_bend_angle.png')
    PLT.close(fig2)

def PlotVelocity(run, outdir):
    """Plot velocity vs time for an individual run.
    
    Parameters:
        run (CalTorque): single CalTorque object to plot
        outdir (str): path to folder where plot will be dropped
    """
    
    vel = run.Velocity()
    
    # get start and end of run
    start = 0
    while (run.NewPosition()[start] == run.NewPosition()[0]):
        start += 1
    end = len(run.NewPosition()) - 1
    while (run.NewPosition()[end] == run.NewPosition()[-1]):
        end -= 1
    
    # slice off flat ends from data
    CHOP_START = start/6
    CHOP_END = end + (len(run.NewPosition()) - end) / 6
    timeChop = run.Time()[CHOP_START : CHOP_END]
    posChop = run.NewPosition()[CHOP_START : CHOP_END]
    velChop = vel[CHOP_START : CHOP_END]
    
    # create plot folder if it does not exist
    if not os.path.exists(outdir):
        os.makedirs(outdir)
    
    # plot
    fig = PLT.figure(figsize = (15, 8), dpi = 150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax = fig.add_axes(box)
    ax2 = ax.twinx()
    ax.set_ylabel('Position (cm)', color = 'blue', fontsize = 18)
    ax.set_xlabel('Time (s)', fontsize = 18)
    ax2.set_ylabel('Velocity (cm/s)', color = 'green', fontsize = 18)
    PLT.title(run.Datetime().strftime('%Y%m%d') + ' Source ' + run.Source() + ' ' + run.Direction() + ' ' + run.Station() + ' Velocity vs Time', fontsize = 20, y = 1.02)
    ax.plot(timeChop, posChop, color = 'blue', linewidth = 2)
    ax2.plot(timeChop, velChop, color = 'green', linewidth = 2)
    ax2.grid()
    ax.tick_params(labelsize=18)
    ax2.tick_params(labelsize=18)
    fig.savefig(outdir + run.Datetime().strftime('%Y%m%d') + run.Filename().split('-20')[-1].split('.')[0][2:] + '_velocity.png')
    PLT.close(fig)
    
def PlotAvgTorque(data, average, dir, bins, max_d, outdir):
    """Plot average torque for all deployments/retractions to/from S5 using weakTh-228
    
    Parameters:
        data (dict): dictionary of new CalTorque objects to add to averages and update plot with
        average (dict): dictionary of previously computed averages, can be empty
        bins (int): number of bins
        max_d (int): maximum displacement undergone by source
        outdir (str): path to folder where plot will be dropped
        
    Returns:
        average updated with the new data
    """
    
    x_inputs = numpy.zeros(bins)
    X = numpy.linspace(0,max_d,bins)
    # if dictionary is empty, create contents
    if any(average) == False:
        average['torAvg'] = numpy.zeros(bins)
        average['stanDev'] = numpy.zeros(bins)
        average['allTorques'] = []
        for i in range(X.size):
            average['allTorques'].append([])
    
    # bin and add new data to average
    for run in data:
       if data[run].Direction() != dir: continue #Skip all retractions
       if data[run].Source() != 'weakTh-228': continue #Skip all other sources
       if 'S5' not in data[run].Station(): continue # Skip calibration runs
       binx, biny = Bin(data[run].NewPosition(), data[run].Torque(), dl = 1, xmin = 0, xmax = 456)
       
       for j, b in enumerate(biny):
          if b == 0: continue
          average['torAvg'][j] += b
          average['allTorques'][j].append(b)
          x_inputs[j]+=1.
    
    # calculate average for each bin
    for i,n in enumerate(x_inputs):
       if n==0: continue
       average['torAvg'][i] = average['torAvg'][i]/n
       average['stanDev'][i] = numpy.std(average['allTorques'][i], dtype = numpy.float64)
    
    # plot average
    fig = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax1 = fig.add_axes(box)
    ax1.set_xlabel('Position (cm)', fontsize = 18)
    ax1.set_ylabel('Torque (in.oz)', fontsize = 18)
    miny = min(average['stanDev'])
    maxy = max(average['torAvg'])
    maxx = max(X)
    
    ## Add station locations
    if maxx >= 328.9:
        ax1.axvline(338.9, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S2', xy = (340.9, miny * 0.8), color = 'r')
    if maxx >= 382.6:
        ax1.axvline(392.6, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S5', xy = (394.6, miny * 0.8), color = 'r')
    if maxx >= 435.7:
        ax1.axvline(445.7, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S8', xy = (447.7, miny * 0.8), color = 'r')
    if maxx >= 491.7:
        ax1.axvline(501.7, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S11', xy = (503.7, miny * 0.8), color = 'r')
    if maxx >= 595.8:
        ax1.axvline(605.8, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S17', xy = (607.8, miny * 0.8), color = 'r')
    
    ## add bend locations
    if maxx >= 47.8:
        ax1.axvspan(47.8, 93.5, color = 'g', alpha = 0.2)
        ax1.annotate('Bellows', xy = (48, miny), color = 'g')
    if maxx >= 136.4:
        ax1.axvspan(136.4, 144.6, color = 'g', alpha = 0.2)
        ax1.annotate('First outer bend', xy = (136, miny), color = 'g')
    if maxx >= 196.5:
        ax1.axvspan(196.5, 218.0, color = 'g', alpha = 0.2)
        ax1.annotate('The "S"', xy = (212, miny), color = 'g')
    if maxx >= 228.3:
        ax1.axvspan(228.3, 241.5, color = 'g', alpha = 0.2)
    if maxx >= 274.6:
        ax1.axvspan(274.6, 284.6, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A5', xy = (274, miny), color = 'g')
    if maxx >= 309.0:
        ax1.axvspan(309.0, 318.9, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A6', xy = (309, miny), color = 'g')
    if maxx >= 358.8:
        ax1.axvspan(358.8, 370.8, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A7', xy = (359, miny), color = 'g')
    if maxx >= 414.5:
        ax1.axvspan(414.5, 426.5, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A8', xy = (414, miny), color = 'g')
    if maxx >= 445.3:
        ax1.axvspan(445.3, 455.3, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A9', xy = (445, miny), color = 'g')
    if maxx >= 467.5:
        ax1.axvspan(467.5, 479.5, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A10', xy = (467, miny), color = 'g')
    if maxx >= 523.8:
        ax1.axvspan(523.8, 535.8, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A11', xy = (523, miny), color = 'g')
    if maxx >= 571.6:
        ax1.axvspan(571.6, 583.6, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A12', xy = (571, miny), color = 'g')
    if maxx >= 627.9:
        ax1.axvspan(627.9, 639.9, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A13', xy = (627, miny), color = 'g')
    if maxx >= 652.7:
        ax1.axvspan(652.7, 660.7, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A14', xy = (652, miny), color = 'g')
    if maxx >= 678.5:
        ax1.axvspan(678.5, 681.1, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A15', xy = (678, miny), color = 'g')
    
    PLT.title("Average of all weak Th-228" + dir + " deployments to S5", fontsize = 20, y = 1.02)
    ax1.grid()
    ax1.tick_params(labelsize=18)
    PLT.errorbar(X, average['torAvg'], xerr = None, yerr = average['stanDev'], linewidth = 2, elinewidth = 1, capsize=0, ecolor = (.4, .4, 1))
    fig.savefig(outdir + "AverageTorque_" + dir + ".png")
    PLT.close(fig)
    return average

def PlotAvgTorHist(avg, dir, outdir):
    """Plot distribution of torques for each point on the average torque deployments plot
    
    Parameters:
        avg (list): average torque data
        dir (str): direction of run (in or out)
        outdir (str): path to folder where plots will be dropped
    """
    for i, a in enumerate(avg):
        lower = int(round(min(a))) - 1
        higher = int(round(max(a))) + 1
        bins = arange(lower, higher, 0.5)
        hist, binEdges = histogram(a, bins)
        fig = PLT.figure(figsize=(15,8),dpi=150)
        box = [0.14, 0.14, 0.76, 0.76]
        ax1 = fig.add_axes(box)
        ax1.set_xlabel('Torque (in.oz)', fontsize = 18)
        ax1.set_ylabel('Number of Runs', fontsize = 18)
        ax1.set_xlim(lower, higher)
        ax1.grid()
        ax1.tick_params(labelsize=18)
        PLT.title("Torque distribution for bin %d of average torque %s"%(i, dir), fontsize = 20, y = 1.02)
        PLT.bar(binEdges[:-1], hist, width = binEdges[1] - binEdges[0])
        PLT.savefig("Histograms/%s/AvgTor_Bin%d.png"%(i,dir))
        PLT.close(fig)

def PlotSources(data, src, bins, max_d, outdir):
    """Plot average torque for each source type.
        
    Parameters:
        data (dict): dictionary of new CalTorque objects to add to averages and update plot with
        month (dict): dictionary of previously computed averages, can be an empty sequence
        bins (int): number of bins
        max_d (int): maximum displacement undergone by source
        outdir (str): path to folder where plot will be dropped
    Returns:
        computed averages 
    """
    
    # store position and torque data for each run in appropriate src entry
    for run in data:
        if data[run].Datetime().strftime('%Y%m%d') <= 20130222: continue
        if data[run].Source() not in src.keys():
            src[data[run].Source()] = {'X': numpy.linspace(0,max_d,bins), 'Y': numpy.zeros(bins), 'numRuns': 0}
    
    # bin new data
    for s in src:
        x_inputs = numpy.zeros(bins)
        X = src[s]['X']
        Y = src[s]['Y']
        counts = src[s]['numRuns']
        for run in data:
            if data[run].Direction() != 'in': continue #Skip all retractions
            if data[run].Datetime().strftime('%Y%m%d') <= 20130222: continue
            if data[run].Source() != s: continue
            binx, biny = Bin(data[run].NewPosition(), data[run].Torque(), max_d/bins)
            counts += 1
            # add binned run data to average data
            for j, b in enumerate(biny):
                if b == 0: continue
                Y[j] += b
                x_inputs[j]+=1.
        # calculate average torque in a bin
        for l, n in enumerate(x_inputs):
            if n==0: continue
            Y[l] = Y[l]/n
        src[s]['X'] = X
        src[s]['Y'] = Y
        src[s]['numRuns'] = counts
    
    # plot source averages
    fig = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax = fig.add_axes(box)
    ax.set_xlabel('Position (cm)', fontsize = 18)
    ax.set_ylabel('Torque (in.oz)', fontsize = 18)
    maxx = max(src['weakCs-137']['X'])
    miny = 0
    
    ## Add station locations
    if maxx >= 328.9:
        ax.axvline(338.9, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S2', xy = (340.9, miny * 0.8), color = 'r')
    if maxx >= 382.6:
        ax.axvline(392.6, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S5', xy = (394.6, miny * 0.8), color = 'r')
    if maxx >= 435.7:
        ax.axvline(445.7, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S8', xy = (447.7, miny * 0.8), color = 'r')
    if maxx >= 491.7:
        ax.axvline(501.7, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S11', xy = (503.7, miny * 0.8), color = 'r')
    if maxx >= 595.8:
        ax.axvline(605.8, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S17', xy = (607.8, miny * 0.8), color = 'r')
    
    ## Add bend locations
    if maxx >= 47.8:
        ax.axvspan(47.8, 93.5, color = 'g', alpha = 0.2)
        ax.annotate('Bellows', xy = (48, miny), color = 'g')
    if maxx >= 136.4:
        ax.axvspan(136.4, 144.6, color = 'g', alpha = 0.2)
        ax.annotate('First outer bend', xy = (136, miny), color = 'g')
    if maxx >= 196.5:
        ax.axvspan(196.5, 218.0, color = 'g', alpha = 0.2)
        ax.annotate('The "S"', xy = (212, miny), color = 'g')
    if maxx >= 228.3:
        ax.axvspan(228.3, 241.5, color = 'g', alpha = 0.2)
    if maxx >= 274.6:
        ax.axvspan(274.6, 284.6, color = 'g', alpha = 0.2)
        ax.annotate('Bend A5', xy = (274, miny), color = 'g')
    if maxx >= 309.0:
        ax.axvspan(309.0, 318.9, color = 'g', alpha = 0.2)
        ax.annotate('Bend A6', xy = (309, miny), color = 'g')
    if maxx >= 358.8:
        ax.axvspan(358.8, 370.8, color = 'g', alpha = 0.2)
        ax.annotate('Bend A7', xy = (359, miny), color = 'g')
    if maxx >= 414.5:
        ax.axvspan(414.5, 426.5, color = 'g', alpha = 0.2)
        ax.annotate('Bend A8', xy = (414, miny), color = 'g')
    if maxx >= 445.3:
        ax.axvspan(445.3, 455.3, color = 'g', alpha = 0.2)
        ax.annotate('Bend A9', xy = (445, miny), color = 'g')
    if maxx >= 467.5:
        ax.axvspan(467.5, 479.5, color = 'g', alpha = 0.2)
        ax.annotate('Bend A10', xy = (467, miny), color = 'g')
    if maxx >= 523.8:
        ax.axvspan(523.8, 535.8, color = 'g', alpha = 0.2)
        ax.annotate('Bend A11', xy = (523, miny), color = 'g')
    if maxx >= 571.6:
        ax.axvspan(571.6, 583.6, color = 'g', alpha = 0.2)
        ax.annotate('Bend A12', xy = (571, miny), color = 'g')
    if maxx >= 627.9:
        ax.axvspan(627.9, 639.9, color = 'g', alpha = 0.2)
        ax.annotate('Bend A13', xy = (627, miny), color = 'g')
    if maxx >= 652.7:
        ax.axvspan(652.7, 660.7, color = 'g', alpha = 0.2)
        ax.annotate('Bend A14', xy = (652, miny), color = 'g')
    if maxx >= 678.5:
        ax.axvspan(678.5, 681.1, color = 'g', alpha = 0.2)
        ax.annotate('Bend A15', xy = (678, miny), color = 'g')
    
    colors = "bgrcmy"
    color_index = 0
    for s in src:
        if src[s]['Y'] == []: continue
        if s == 'unknown': continue
        # get rid of empty bins before plotting
        noZero = numpy.nonzero(src[s]['Y'])# return tuple of an array w/ indices of nonZero values
        Y_nonull = []
        X_nonull = []
        for i, x in enumerate(src[s]['X']):
            if i in noZero[0]:
                Y_nonull.append(src[s]['Y'][i])
                X_nonull.append(src[s]['X'][i])
        # plot source
        if s == 'weakTh-228':
            ax.plot(X_nonull, Y_nonull, '-', linewidth = 2, c='k', label = "weakTh-228, " + str(src[s]['numRuns']) + " runs")
        else:
            ax.plot(X_nonull, Y_nonull, '-', linewidth = 2, c=colors[color_index], label = s + ", " + str(src[s]['numRuns']) + " runs")
            color_index += 1
    h, l = ax.get_legend_handles_labels()
    ax.grid()
    ax.tick_params(labelsize=18)
    PLT.title('Average Torque for Different Sources', fontsize = 20, y = 1.02)
    PLT.legend(h, l, 'upper left')
    PLT.savefig(outdir + 'AvgTorqueSources.png')
    PLT.close(fig)
    return src


def PlotMonthlyAvgs(data, months, bins, max_d, outdir):
    """Plot average torque for all deployments in a month using weakTh-228.
    
    Parameters:
        data (dict): dictionary of new CalTorque objects to add to averages and update plot with
        month (dict): dictionary of previously computed averages, can be an empty sequence
        bins (int): number of bins
        max_d (int): maximum displacement undergone by source
        outdir (str): path to folder where plot will be dropped
        
    Returns:
        monthly averages updated with the new data
    """
    
    # retrieve new month entries and store as keys
    for run in data:
        month_year = data[run].Datetime().strftime('%b%Y')
        if month_year not in months.keys():
            months[month_year] = {'X': numpy.linspace(0,max_d,bins), 'Y': numpy.zeros(bins), 'Y2': []}
            for i in range(months[month_year]['X'].size):
                months[month_year]['Y2'].append([])
    
    # bin and average data
    for m in months:
        x_inputs = numpy.zeros(bins)
        X = months[m]['X']
        Y = months[m]['Y']
        Y2 = months[m]['Y2']
        # append data of each file to appropriate month dictionary entry
        for run in data:
            if data[run].Direction() != 'in': continue #Skip all retractions
            if data[run].Datetime().strftime('%b%Y') not in m: continue
            if data[run].Source() != 'weakTh-228': continue
            if data[run].Station() in ['S8', 'S11', 'Other']: continue
            binx, biny = Bin(data[run].NewPosition(), data[run].Torque(), max_d/bins)
            
            for j, b in enumerate(biny):
                if b == 0: continue
                Y[j] += b
                Y2[j].append(b)
                x_inputs[j]+=1.
        Y_std = numpy.zeros(len(X))
        # average content of each bin, calculate standard deviation
        for l, n in enumerate(x_inputs):
            if n!=0: Y[l] = Y[l]/n
            Y_std[l] = numpy.std(Y2[l])
        months[m]['X'] = X
        months[m]['Y'] = Y
        months[m]['Y2'] = Y2
        
        #_____plot individual monthly averages)_____#
        
        # checks in folder where plots are dropped to see if data has already been plotted
        alreadyPlotted = []
        dropFolder = os.listdir("averagePlots")
        for file in dropFolder:
            alreadyPlotted.append(file.split('.')[0].split('_')[-1])
        for run in data:
            if data[run].Datetime().strftime('%m%Y') in alreadyPlotted:
                continue

        # plot
        fig = PLT.figure(figsize=(15,8),dpi=150)
        box = [0.14, 0.14, 0.76, 0.76]
        ax1 = fig.add_axes(box)
        ax1.set_xlabel('Position (cm)', fontsize = 18)
        ax1.set_ylabel('Torque (in.oz)', fontsize = 18)
        miny = min(Y_std)
        maxy = max(Y)
        maxx = max(X)
        
        ## Add station locations
        if maxx >= 328.9:
            ax1.axvline(338.9, color = 'r', linewidth = 1, linestyle = '-')
            ax1.annotate('S2', xy = (340.9, miny * 0.8), color = 'r')
        if maxx >= 382.6:
            ax1.axvline(392.6, color = 'r', linewidth = 1, linestyle = '-')
            ax1.annotate('S5', xy = (394.6, miny * 0.8), color = 'r')
        if maxx >= 435.7:
            ax1.axvline(445.7, color = 'r', linewidth = 1, linestyle = '-')
            ax1.annotate('S8', xy = (447.7, miny * 0.8), color = 'r')
        if maxx >= 491.7:
            ax1.axvline(501.7, color = 'r', linewidth = 1, linestyle = '-')
            ax1.annotate('S11', xy = (503.7, miny * 0.8), color = 'r')
        if maxx >= 595.8:
            ax1.axvline(605.8, color = 'r', linewidth = 1, linestyle = '-')
            ax1.annotate('S17', xy = (607.8, miny * 0.8), color = 'r')
        
        ## Add bend locations
        if maxx >= 47.8:
            ax1.axvspan(47.8, 93.5, color = 'g', alpha = 0.2)
            ax1.annotate('Bellows', xy = (48, miny), color = 'g')
        if maxx >= 136.4:
            ax1.axvspan(136.4, 144.6, color = 'g', alpha = 0.2)
            ax1.annotate('First outer bend', xy = (136, miny), color = 'g')
        if maxx >= 196.5:
            ax1.axvspan(196.5, 218.0, color = 'g', alpha = 0.2)
            ax1.annotate('The "S"', xy = (212, miny), color = 'g')
        if maxx >= 228.3:
            ax1.axvspan(228.3, 241.5, color = 'g', alpha = 0.2)
        if maxx >= 274.6:
            ax1.axvspan(274.6, 284.6, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A5', xy = (274, miny), color = 'g')
        if maxx >= 309.0:
            ax1.axvspan(309.0, 318.9, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A6', xy = (309, miny), color = 'g')
        if maxx >= 358.8:
            ax1.axvspan(358.8, 370.8, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A7', xy = (359, miny), color = 'g')
        if maxx >= 414.5:
            ax1.axvspan(414.5, 426.5, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A8', xy = (414, miny), color = 'g')
        if maxx >= 445.3:
            ax1.axvspan(445.3, 455.3, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A9', xy = (445, miny), color = 'g')
        if maxx >= 467.5:
            ax1.axvspan(467.5, 479.5, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A10', xy = (467, miny), color = 'g')
        if maxx >= 523.8:
            ax1.axvspan(523.8, 535.8, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A11', xy = (523, miny), color = 'g')
        if maxx >= 571.6:
            ax1.axvspan(571.6, 583.6, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A12', xy = (571, miny), color = 'g')
        if maxx >= 627.9:
            ax1.axvspan(627.9, 639.9, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A13', xy = (627, miny), color = 'g')
        if maxx >= 652.7:
            ax1.axvspan(652.7, 660.7, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A14', xy = (652, miny), color = 'g')
        if maxx >= 678.5:
            ax1.axvspan(678.5, 681.1, color = 'g', alpha = 0.2)
            ax1.annotate('Bend A15', xy = (678, miny), color = 'g')
        
        PLT.title('Average of '+ m + ' IN deployments', fontsize = 20, y = 1.02)
        ax1.grid()
        ax1.tick_params(labelsize=18)
        PLT.errorbar(X, Y, xerr = None, yerr = Y_std, linewidth = 2, elinewidth = 1, capsize=0, ecolor = (.4, .4, 1))
        fig.savefig(outdir + 'average_in_' + m + '.png')
    
    #_____plot all monthly averages together_____#
    fig = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax = fig.add_axes(box)
    ax.set_xlabel('Position (cm)', fontsize = 18)
    ax.set_ylabel('Torque (in.oz)', fontsize = 18)
    colorthis = [0, 0, 1]
    denom = int(len(months)) - 6.
    for i, month in enumerate(sorted(months, key = lambda x:datetime.strptime(x,'%b%Y'))):
        if month == 'Oct2012' or month == 'Dec2012' or month == 'Nov2012' or month == 'Jan2013' or month == 'Feb2013' or month == 'Mar2012': continue #  or month == 'Oct2013' because October only has one run
        
        # get rid of empty bins before plotting
        noZero = numpy.nonzero(months[m]['Y'])# return tuple of an array w/ indices of nonZero values
        Y_nonull = []
        X_nonull = []
        Y2_nonull = []
        for i, x in enumerate(months[m]['X']):
            if i in noZero[0]:
                Y_nonull.append(months[m]['Y'][i])
                X_nonull.append(months[m]['X'][i])
                Y2_nonull.append(months[m]['Y2'][i])
        ax.plot(months[month]['X'], months[month]['Y'], '-', linewidth = 1, c = tuple(colorthis), label = month)
        colorthis[0] += (1 / denom)
        colorthis[2] -= (1 / denom)
    h, l = ax.get_legend_handles_labels()
    ax.grid()
    ax.tick_params(labelsize=18)
    PLT.title('Monthly Torque Average (after Feb 2013)', fontsize = 20, y = 1.02)
    PLT.legend(h, l, 'upper left', ncol = 2)
    PLT.savefig(outdir + 'Monthly_averages.png')
    PLT.close(fig)
    return months
    
def PlotTelescope(data, telescope, bins, max_d, outdir):
    """Averages torque for runs before or after Feb 22, 2013. Plots both averages together
        for comparison.
            
    Parameters:
        data (dict): dictionary of new CalTorque objects to add to averages and update plot with
        telescope (dict): dictionary of previously computed averages, can be an empty sequence
        bins (int): number of bins
        max_d (int): maximum displacement undergone by source
        outdir (str): path to folder where plot will be dropped
    
    Returns:
        telescope, average updated with the new data
    """
    
    # create telescope data structure if does not exist
    if any(telescope) == False:
        telescope['before'] = [numpy.linspace(0,max_d,bins), numpy.zeros(bins)]
        telescope['after'] = [numpy.linspace(0,max_d,bins), numpy.zeros(bins)]
    
    # bin and average data
    for t in telescope:
        x_inputs = numpy.zeros(bins)
        X = telescope[t][0]
        Y = telescope[t][1]
        for run in data:
            if data[run].Direction() == 'out': continue
            if data[run].Source() != 'weakTh-228': continue
            if data[run].Station() in ['S8', 'S11', 'Other']: continue
            date = int(data[run].Datetime().strftime('%Y%m%d'))
            if date <= 20130222 and t != 'before': continue
            elif date > 20130222 and t != 'after': continue
            binx, biny = Bin(data[run].NewPosition(), data[run].Torque(), max_d/bins)
            # add data to average
            for j, b in enumerate(biny):
                if b == 0.: continue
                Y[j] += b
                x_inputs[j]+=1.
        # calculate average
        for l, n in enumerate(x_inputs):
            if n == 0: continue
            Y[l] = Y[l]/n
        telescope[t] = [X, Y]
        
    # make plot
    fig = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax = fig.add_axes(box)
    ax.set_xlabel('Position (cm)', fontsize = 18)
    ax.set_ylabel('Torque (in.oz)', fontsize = 18)
    
    miny = 1
    maxx = 400
    
    ## add box at telescope location
    ax.axvspan(20.8, 33.3, color = 'magenta', alpha = 0.2)
    ax.annotate('Telescope Travel', xy = (20.8, miny + 2), color = 'magenta')
            
    ## Add station locations
    if maxx >= 328.9:
        ax.axvline(338.9, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S2', xy = (340.9, miny * 0.8), color = 'r')
    if maxx >= 382.6:
        ax.axvline(392.6, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S5', xy = (394.6, miny * 0.8), color = 'r')
    if maxx >= 435.7:
        ax.axvline(445.7, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S8', xy = (447.7, miny * 0.8), color = 'r')
    if maxx >= 491.7:
        ax.axvline(501.7, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S11', xy = (503.7, miny * 0.8), color = 'r')
    if maxx >= 595.8:
        ax.axvline(605.8, color = 'r', linewidth = 1, linestyle = '-')
        ax.annotate('S17', xy = (607.8, miny * 0.8), color = 'r')
    
    ## Add bend locations
    if maxx >= 47.8:
        ax.axvspan(47.8, 93.5, color = 'g', alpha = 0.2)
        ax.annotate('Bellows', xy = (48, miny), color = 'g')
    if maxx >= 136.4:
        ax.axvspan(136.4, 144.6, color = 'g', alpha = 0.2)
        ax.annotate('First outer bend', xy = (136, miny), color = 'g')
    if maxx >= 196.5:
        ax.axvspan(196.5, 218.0, color = 'g', alpha = 0.2)
        ax.annotate('The "S"', xy = (212, miny), color = 'g')
    if maxx >= 228.3:
        ax.axvspan(228.3, 241.5, color = 'g', alpha = 0.2)
    if maxx >= 274.6:
        ax.axvspan(274.6, 284.6, color = 'g', alpha = 0.2)
        ax.annotate('Bend A5', xy = (274, miny), color = 'g')
    if maxx >= 309.0:
        ax.axvspan(309.0, 318.9, color = 'g', alpha = 0.2)
        ax.annotate('Bend A6', xy = (309, miny), color = 'g')
    if maxx >= 358.8:
        ax.axvspan(358.8, 370.8, color = 'g', alpha = 0.2)
        ax.annotate('Bend A7', xy = (359, miny), color = 'g')
    if maxx >= 414.5:
        ax.axvspan(414.5, 426.5, color = 'g', alpha = 0.2)
        ax.annotate('Bend A8', xy = (414, miny), color = 'g')
    if maxx >= 445.3:
        ax.axvspan(445.3, 455.3, color = 'g', alpha = 0.2)
        ax.annotate('Bend A9', xy = (445, miny), color = 'g')
    if maxx >= 467.5:
        ax.axvspan(467.5, 479.5, color = 'g', alpha = 0.2)
        ax.annotate('Bend A10', xy = (467, miny), color = 'g')
    if maxx >= 523.8:
        ax.axvspan(523.8, 535.8, color = 'g', alpha = 0.2)
        ax.annotate('Bend A11', xy = (523, miny), color = 'g')
    if maxx >= 571.6:
        ax.axvspan(571.6, 583.6, color = 'g', alpha = 0.2)
        ax.annotate('Bend A12', xy = (571, miny), color = 'g')
    if maxx >= 627.9:
        ax.axvspan(627.9, 639.9, color = 'g', alpha = 0.2)
        ax.annotate('Bend A13', xy = (627, miny), color = 'g')
    if maxx >= 652.7:
        ax.axvspan(652.7, 660.7, color = 'g', alpha = 0.2)
        ax.annotate('Bend A14', xy = (652, miny), color = 'g')
    if maxx >= 678.5:
        ax.axvspan(678.5, 681.1, color = 'g', alpha = 0.2)
        ax.annotate('Bend A15', xy = (678, miny), color = 'g')
    
    ax.plot(telescope['before'][0], telescope['before'][1], '-', linewidth = 2, c = 'b', label = 'Before')
    ax.plot(telescope['after'][0], telescope['after'][1], '-', linewidth = 2, c = 'r', label = 'After')
    h, l = ax.get_legend_handles_labels()
    ax.grid()
    ax.tick_params(labelsize=18)
    PLT.title('Average Torque Before and After Feb 22, 2013', fontsize = 20, y = 1.02)
    PLT.legend(h, l, 'upper left')
    PLT.savefig(outdir + 'Before_After_New_Telescope.png')
    
    return telescope

# _________ Plot Max torque vs date _________#
def Collect(run):
    """Retrieve position and torque data from database entry.
    
    Parameters:
        run (CalTorque): single CalTorque object
    
    Returns:
        maximum torque for the entire run, for the area around Bend A6, and for the area around Bend A7
    """
    torque = run.Torque()
    position = run.Position()
    maxTorque = max(torque)
    
    pos_torque = []
    for i, x in enumerate(position):
        if x > 300 and x < 330:
            pos_torque.append(torque[i])
    max_at_B2 = (max(pos_torque))
    
    pos_torque2 = []
    for i, x in enumerate(position):
        if x > 350 and x < 375:
            pos_torque2.append(torque[i])
    max_at_B3 = (max(pos_torque2))
    
    return maxTorque, max_at_B2, max_at_B3

def Linearfnc(x, m, b):
    """Linear function model
    
    Parameters:
        x (float): position variable
        m (float): slope
        b (float): coefficient
    """
    return m*x + b
    
def LinearFit(date, maxTorque):
    """Perform linear fit to max torque vs date data 
    
    Parameters:
        date (int, array): array of date values from each run
        maxTorque (float, array): array of maximum torque values from each run
        
    Returns:
        fit parameters (slope and coefficient)
    """
    fit = linregress(date, maxTorque)
    return fit


def FitArray(fit, date):
    """Produce array of y values from linear fit equation parameters
    
    Parameters:
        fit (float, array): linear fit parameters (slope and coefficient)
        date (int, array): array of date values from each run
        
    Returns:
        y values of linear fit function
    """
    fitArray = []
    for day in date:
        fitArray.append(Linearfnc(day, fit[0], fit[1]))
    return fitArray


def PlotMaxtvsDate(data, outdir):
    """Plot maximum torque vs date for each run
    
    Parameters:
        data (dict): dictionary of CalTorque objects
        outdir (str): path to folder where plot will be dropped
    """
    max_torque = []
    max_torque_Before = []
    max_at_B2 = []
    max_at_B2_Before = []
    max_at_B3 = []
    max_at_B3_Before = []
    date = []
    date_Before = []
    intdate = []
    intdate_Before = []
    filelist = []
    filesBefore = []
    files = []
    linfitBefore = []
    linfit = []
    linfit_B2_Before = []
    linfit_B3_Before = []
    linfit_B2 = []
    linfit_B3 = []
    
    for run in data:
        if data[run].Direction() != 'in': continue
        if data[run].Source() != 'weakTh-228': continue
        if 'S5' not in data[run].Station(): continue
        if max(data[run].NewPosition()) < 350: continue
        if data[run].NewPosition()[0] > 50: continue
        if data[run].Filename() == 'Live Graph 20-Mar-2012_?.xlsx': continue
        try:
            maxtorque, maxB2, maxB3 = Collect(data[run])
        except:
            print data[run].Filename(), "Sequence was empty. Skipped run."
            continue
        intd = int(data[run].Datetime().strftime('%Y%m%d'))
        
        # sort torque values based on date of deployment
        if intd <= 20130222:
            max_at_B2_Before.append(maxB2)
            max_at_B3_Before.append(maxB3)
            max_torque_Before.append(maxtorque)
            date_Before.append(data[run].Datetime())
        else:
            max_at_B2.append(maxB2)
            max_at_B3.append(maxB3)
            max_torque.append(maxtorque)
            date.append(data[run].Datetime())
            
    # convert dates to ints for construction of plots
    for d in date:
        intd = d - min(date)
        intdate.append(intd.days)
    for d in date_Before:
        intd = d - min(date_Before)
        intdate_Before.append(intd.days)
    
    # fit linear trendlines to torque vs date points
    lin = LinearFit(intdate, max_torque)
    lin_B2 = LinearFit(intdate, max_at_B2)
    lin_B3 = LinearFit(intdate, max_at_B3)
    lin_B = LinearFit(intdate_Before, max_torque_Before)
    lin_B2_B = LinearFit(intdate_Before, max_at_B2_Before)
    lin_B3_B = LinearFit(intdate_Before, max_at_B3_Before)
    
    # get y values for the fit
    linfit = FitArray(lin, intdate)
    linfit_B2 = FitArray(lin_B2, intdate)
    linfit_B3 = FitArray(lin_B3, intdate)
    linfit_B = FitArray(lin_B, intdate_Before)
    linfit_B2_B = FitArray(lin_B2_B, intdate_Before)
    linfit_B3_B = FitArray(lin_B3_B, intdate_Before)
    
    # plot
    fig = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax1 = fig.add_axes(box)
    ax1.set_ylabel('Max Torque (in.oz)', fontsize = 18)
    ax1.set_xlabel('Date', fontsize = 18)
    PLT.title('Max Torque Versus Date', fontsize = 20, y = 1.02)
    ax1.grid()
    ax1.tick_params(labelsize=18)
    ax1.scatter(date, max_torque, color = 'g', label = 'Data')
    ax1.plot(date, linfit, '-', color = 'g', label = 'Linear fit: y = %.2f * x + %.2f' % (lin[0], lin[1]))
    ax1.scatter(date_Before, max_torque_Before, color = 'g')
    ax1.plot(date_Before, linfit_B, '-', color = 'black', label = 'Linear fit: y = %.2f * x + %.2f' % (lin_B[0], lin_B[1]))
    #ax1.axvline(datetime.strptime("2013-02-28", "%Y-%m-%d"), colors = 'magenta')
    ax1.text(15, 8, "B refers to before the telescope break", size = 'small')
    h, l = ax1.get_legend_handles_labels()
    PLT.legend(h, l, 'lower left')
    fig.savefig(outdir + 'MaxTorque_vs_Date.png')
    PLT.close(fig)

    fig2 = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax2 = fig2.add_axes(box)
    ax2.set_ylabel('Max Torque (in.oz)', fontsize = 18)
    ax2.set_xlabel('Date', fontsize = 18)
    PLT.title('Max Torque Versus Date for the Second Inner Bend', fontsize = 20, y = 1.02)
    ax2.grid()
    ax2.tick_params(labelsize=18)
    ax2.scatter(date, max_at_B2, color = 'b', label = 'Data')
    ax2.plot(date, linfit_B2, '-', color = 'b', label = 'Linear fit: y = %.2f * x + %.2f' % (lin_B2[0], lin_B2[1]))
    ax2.scatter(date_Before, max_at_B2_Before, color = 'b')
    ax2.plot(date_Before, linfit_B2_B, '-', color = 'black', label = 'Linear fit: y = %.2f * x + %.2f' % (lin_B2_B[0], lin_B2_B[1]))
    ax2.annotate("B refers to before the telescope break", xy=(0.9, 0), size = 'small')
    h, l = ax2.get_legend_handles_labels()
    PLT.legend(h, l, 'lower left')
    fig2.savefig(outdir + 'MaxTorque_vs_Date_for_B2.png')
    PLT.close(fig2)
    
    fig3 = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax3 = fig3.add_axes(box)
    ax3.set_ylabel('Max Torque (in.oz)', fontsize = 18)
    ax3.set_xlabel('Date', fontsize = 18)
    PLT.title('Max Torque Versus Date for the Third Inner Bend', fontsize = 20, y = 1.02)
    ax3.grid()
    ax3.tick_params(labelsize=18)
    ax3.scatter(date, max_at_B3, color = 'r', marker = 'D', label = 'Data')
    ax3.plot(date, linfit_B3, '-', color = 'r', label = 'Linear fit: y = %.2f * x + %.2f' % (lin_B3[0], lin_B3[1]))
    ax3.scatter(date_Before, max_at_B3_Before, color = 'r', marker = 'D')
    ax3.plot(date_Before, linfit_B3_B, '-', color = 'black', label = 'Linear fit: y = %.2f * x + %.2f' % (lin_B3_B[0], lin_B3_B[1]))
    ax3.annotate("B refers to before the telescope break", xy=(0.9, 0), size = 'small')
    h, l = ax3.get_legend_handles_labels()
    PLT.legend(h, l, 'lower left')
    fig3.savefig(outdir + 'MaxTorque_vs_Date_for_B3.png')
    PLT.close(fig3)
    
    fig4 = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax4 = fig4.add_axes(box)
    ax4.set_ylabel('Max Torque (in.oz)', fontsize = 18)
    ax4.set_xlabel('Date', fontsize = 18)
    PLT.title('Max Torque Versus Date for Bends A6 and A7', fontsize = 20, y = 1.02)
    ax4.grid()
    ax4.tick_params(labelsize=18)
    ax4.scatter(date, max_at_B3, color = 'r', marker = 'D', label = 'Bend A7')
    ax4.plot(date, linfit_B3, '-', color = 'r', label = 'Bend A7 linear fit: y = %.2f * x + %.2f' % (lin_B3[0], lin_B3[1]))
    ax4.scatter(date_Before, max_at_B3_Before, color = 'r', marker = 'D')
    ax4.plot(date_Before, linfit_B3_B, '-', color = 'black', label = 'Bend A7 linear fit: y = %.2f * x + %.2f' % (lin_B3_B[0], lin_B3_B[1]))
    ax4.scatter(date, max_at_B2, color = 'b', label = 'Bend A6')
    ax4.plot(date, linfit_B2, '-', color = 'b', label = 'Bend A6 linear fit: y = %.2f * x + %.2f' % (lin_B2[0], lin_B2[1]))
    ax4.scatter(date_Before, max_at_B2_Before, color = 'b')
    ax4.plot(date_Before, linfit_B2_B, '-', color = 'g', label = 'Bend A6 linear fit: y = %.2f * x + %.2f' % (lin_B2_B[0], lin_B2_B[1]))
    #ax4.axvline(datetime.strptime("2013-02-28", "%Y-%m-%d"), colors = 'magenta')
    PLT.text(0.5, 0.5, "B refers to before the telescope break", size = 'small')
    h, l = ax4.get_legend_handles_labels()
    PLT.legend(h, l, 'lower left')
    fig4.savefig(outdir + 'MaxTorque_vs_Date_for_peaks.png')
    PLT.close(fig4)

def PlotMaxtvsVel(data, bins, max_d, outdir):
    """Plot maximum torque vs velocity for all deployments to S5 using weakTh-228.
    
    Parameters:
        data (dict): dictionary of CalTorque objects
        bins (int): number of bins
        max_d (int): maximum displacement undergone by source
        outdir (str): path to folder where plot will be dropped
    """
    max_torque = []
    max_vel = []
    max_at_B2 = []
    max_at_B3 = []
    filelist = []
    avgtor = []
    avgvel = []
    linfit = []
    linfit_B2 = []
    linfit_B3 = []
    linfit2 = []
    linfit3 = []
    
    for run in data:
        if data[run].Direction() != 'in': continue
        if data[run].Source() != 'weakTh-228': continue
        if 'S5' not in data[run].Station(): continue
        if max(data[run].NewPosition()) < 350: continue
        if max(data[run].NewPosition()) > 400: continue
        if data[run].NewPosition()[0] > 50: continue
        
        filelist.append(data[run].Filename().split('.')[0].split(' ')[2].split('_')[0])
        max_torque.append(data[run].MaxTorque())
        max_vel.append(max(data[run].Velocity()))
        
        # average torque
        t = 0
        for torque in data[run].Torque():
            t += torque
        t2 = t / len(data[run].Torque())
        avgtor.append(t2)
        
        # average velocity
        sumvel = 0
        for j, v in enumerate(data[run].Velocity()):
            sumvel += v
        averagevel = sumvel / int(len(data[run].Velocity()))
        avgvel.append(averagevel)
       
        pos_torque = []
        for i, x in enumerate(data[run].NewPosition()):
            if x > 300 and x < 330:
                pos_torque.append(data[run].Torque()[i])
        max_at_B2.append(max(pos_torque))
        
        pos_torque2 = []
        for i, x in enumerate(data[run].NewPosition()):
            if x > 350 and x < 375:
                pos_torque2.append(data[run].Torque()[i])
        max_at_B3.append(max(pos_torque2))
        
    
    lin = linregress(avgvel, max_torque)
    lin_B2 = linregress(avgvel, max_at_B2)
    lin_B3 = linregress(avgvel, max_at_B3)
    lin2 = linregress(avgvel, avgtor)
    lin3 = linregress(max_vel, max_torque)
    
    for thing in avgvel:
        linfit.append(Linearfnc(thing, lin[0], lin[1]))
        linfit_B2.append(Linearfnc(thing, lin_B2[0], lin_B2[1]))
        linfit_B3.append(Linearfnc(thing, lin_B3[0], lin_B3[1]))
        linfit2.append(Linearfnc(thing, lin2[0], lin2[1]))
    for other in max_vel:
        linfit3.append(Linearfnc(other, lin3[0], lin3[1]))
    
    # plot
    box = [0.14, 0.14, 0.76, 0.76]
    
    fig = PLT.figure(figsize=(15,8),dpi=150)
    ax1 = fig.add_axes(box)
    ax1.set_ylabel('Max Torque (in.oz)', fontsize = 18)
    ax1.set_xlabel('Average Velocity (cm/s)', fontsize = 18)
    PLT.title('Max Torque Versus Average Velocity', fontsize = 20, y = 1.02)
    ax1.scatter(avgvel, max_torque, color = 'g')
    ax1.plot(avgvel, linfit, '-', color = 'g')
    ax1.grid()
    ax1.tick_params(labelsize=18)
    #for i, txt in enumerate(filelist):
        #ax1.annotate(txt, (date[i], max_torque[i]), position = (date[i], max_torque[i]), size = 'x-small', alpha = 0.5, rotation = 45)
    fig.savefig(outdir + 'MaxTorque_vs_Velocity.png')
    PLT.close(fig)
    
    fig2 = PLT.figure(figsize=(15,8),dpi=150)
    ax2 = fig2.add_axes(box)
    ax2.set_ylabel('Max Torque (in.oz)', fontsize = 18)
    ax2.set_xlabel('Date', fontsize = 18)
    PLT.title('Max Torque Versus Average Velocity for Bend A6 and A7', fontsize = 20, y = 1.02)
    ax2.grid()
    ax2.tick_params(labelsize=18)
    ax2.scatter(avgvel, max_at_B2, color = 'b', label = 'Bend A6')
    ax2.plot(avgvel, linfit_B2, '-', color = 'b')
    ax2.scatter(avgvel, max_at_B3, color = 'r', marker = 'D', label = 'Bend A7')
    ax2.plot(avgvel, linfit_B3, '-', color = 'r')
    fig2.savefig(outdir + 'MaxTorque_vs_Velocity_for_A67.png')
    PLT.close(fig2)
    
    fig4 = PLT.figure(figsize=(15,8),dpi=150)
    ax4 = fig4.add_axes(box)
    ax4.set_ylabel('Average Torque (in.oz)', fontsize = 18)
    ax4.set_xlabel('Average Velocity (cm/s)', fontsize = 18)
    ax4.grid()
    ax4.tick_params(labelsize=18)
    PLT.title('Average Torque Versus Average Velocity', fontsize = 20, y = 1.02)
    ax4.scatter(avgvel, avgtor, color = 'b')
    ax4.plot(avgvel, linfit2, '-', color = 'b')
    fig4.savefig(outdir + 'AverageTorque_vs_Velocity.png')
    PLT.close(fig4)
    
    fig5 = PLT.figure(figsize=(15,8),dpi=150)
    ax5 = fig5.add_axes(box)
    ax5.set_ylabel('Maximum Torque (in.oz)', fontsize = 18)
    ax5.set_xlabel('Maximum Velocity (cm/s)', fontsize = 18)
    PLT.title('Maximum Torque Versus Maximum Velocity', fontsize = 20, y = 1.02)
    ax5.grid()
    ax5.tick_params(labelsize=18)
    ax5.scatter(max_vel, max_torque, color = 'g')
    fig5.savefig(outdir + 'MaxTorque_vs_MaxVelocity.png')
    PLT.close(fig5)

def PlotTorquevsVelocity(run, outdir):
    """ Plot torque over angular velocity vs angular velocity
    
    Parameters:
        run (dict), data for a single run
        outdir (str), path to folder where plot will be dropped
    """
    tempPos = []
    tempY = []
    tempV = []
    for i, t in enumerate(run.Torque()):
        if run.AngularVelocity()[i] == 0: continue
        y = t / run.AngularVelocity()[i]
        tempY.append(y)
        tempV.append(run.AngularVelocity()[i])
    
    # plot
    box = [0.14, 0.14, 0.76, 0.76]
    fig = PLT.figure(figsize=(15,8),dpi=150)
    ax1 = fig.add_axes(box)
    ax1.set_ylabel("Torque / Angular Velocity (in.oz / rpm)", fontsize = 18)
    ax1.set_xlabel("Angular Velocity (rpm)", fontsize = 18)
    PLT.title("Torque over Angular Velocity vs Angular Velocity", fontsize = 20, y = 1.02)
    ax1.grid()
    ax1.tick_params(labelsize=18)
    ax1.scatter(tempV, tempY, color = [0,0,0.5])
    fig.savefig(outdir + run.Datetime().strftime('%Y%m%d') + run.Filename().split('-20')[-1].split('.')[0][2:] + '_TorOverAngVel_vs_AngVel.png')
    PLT.close(fig)

def PlotOffsetvsP0(data, outdir):
    """Plot offset found using Poly2Gauss vs the initial position for all deployments to S5.
    
    Parameters:
        data (dict): dictionary of CalTorque objects
        outdir (str): path to folder where plot will be dropped
    """
    p0 = []
    offset_G = []

    filelist = []

    for run in data:
        if data[run].Direction() != 'in': continue
        if data[run].Offset() == 0: continue
        if data[run].Station() != 'S5': continue
        if data[run].Position()[0] > 50: continue
        if max(data[run].NewPosition()) < 350: continue
        
        p0.append(data[run].Position()[0])
        filelist.append(data[run].Filename().split('.')[0].split(' ')[2].split('_')[0])
        offset_G.append(data[run].Offset())

    x = list(xrange(-25, 15))
    y = deepcopy(x)

    y1 = []
    y2 = []
    for item in y:
        y1.append(item + 1)
        y2.append(item - 1)
    
    # plot
    fig1 = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax1 = fig1.add_axes(box)
    ax1.set_ylabel('Initial Position Offset (cm)', fontsize = 18)
    ax1.set_xlabel('Poly-Two-Gauss Offset (cm)', fontsize = 18)
    PLT.title('Initial Position Versus Poly-Two-Gauss Offset)', fontsize = 20, y = 1.02)
    ax1.annotate('y = x', xy = (15, 10), color = 'r', fontstyle = 'italic', fontsize = 'x-large')
    ax1.plot(x, y1, '-', color = 'g')
    ax1.plot(x, y2, '-', color = 'g')
    ax1.scatter(offset_G, p0, color = 'maroon')
    #for i, txt in enumerate(filelist): ## NOTE ## Uncomment these lines to annotate data points
        #ax1.annotate(txt, (offset_G[i], p0[i]), textcoords = 'offset points', ha = 'right', va = 'bottom', size = 'x-small', alpha = 0.5, rotation = 45)
    ax1.plot(x, y, '-', color = 'r', linewidth = 2)
    ax1.grid()
    ax1.tick_params(labelsize=18)
    fig1.savefig(outdir + 'p0_vs_Offset_Gauss.png')
    PLT.close(fig1)
    
    thislist = []
    for i, thing in enumerate(p0):
        this = thing - offset_G[i]
        thislist.append(this)
    stan = numpy.std(thislist)

    fig2 = PLT.figure(figsize = (15,8), dpi = 150)
    ax2 = fig2.add_axes(box)
    ax2.set_ylabel('Initial Position spread', fontsize = 18)
    ax2.set_xlabel('Data Offset (cm)', fontsize = 18)
    PLT.title('Initial Position Spread versus Offset (model method)', fontsize = 20, y = 1.02)
    ax2.scatter(offset_G, thislist, color = 'r', label = "Standard Deviation: " + str(stan))
    ax2.grid()
    ax2.tick_params(labelsize=18)
    h, l = ax2.get_legend_handles_labels()
    PLT.legend(h, l, 'lower left')
    fig2.savefig(outdir + 'p0_vs_Offset_Gauss_spread.png')
    PLT.close(fig2)
    
def PlotParameters(data, outdir):
    """ Plot each gaussian fit parameter vs date for all deployments to S5 using weakTh-228.
    
    Fit parameters:
         p[0], p[1], p[2] = amplitude , center_x, sigma: for peak A
         p[3], p[4] = amplitude, sigma: for peak B
         p[5] = polynomial offset
         p[6], p[7], p[8] = ax^2 + bx + c 
    
    Parameters:
        data (dict): dictionary of CalTorque objects
        outdir (str): path to folder where plot will be dropped
    """
    parameters = {}
    parameters['p0'] = []
    parameters['p1'] = []
    parameters['p2'] = []
    parameters['p3'] = []
    parameters['p4'] = []
    parameters['p5'] = []
    parameters['p6'] = []
    parameters['p7'] = []
    parameters['p8'] = []
    filelist = []
    date = []
    intdate = []
        
    # get parameters for each run
    for run in data:
        if data[run].Direction() != 'in': continue
        if data[run].Source() != 'weakTh-228': continue
        if data[run].Station() != 'S5': continue
        if data[run].NewParms() == None: continue
        
        np = data[run].NewParms()
        
        parameters['p0'].append(np[0])
        parameters['p1'].append(np[1])
        parameters['p2'].append(np[2])
        parameters['p3'].append(np[3])
        parameters['p4'].append(np[4])
        parameters['p5'].append(np[5])
        parameters['p6'].append(np[6])
        parameters['p7'].append(np[7])
        parameters['p8'].append(np[8])
        
        filelist.append(data[run].Filename().split('.')[0].split(' ')[2].split('_')[0])
        rundate = datetime.strptime(str(data[run].Datetime()).split(' ')[0], '%Y-%m-%d')
        date.append(rundate)
    
    # convert date to int for plotting
    for d in date:
        intd = d - min(date)
        intdate.append(intd.days)
    
    # fit trendline to each parameter
    l0 = linregress(intdate, parameters['p0'])
    l1 = linregress(intdate, parameters['p1'])
    l2 = linregress(intdate, parameters['p2'])
    l3 = linregress(intdate, parameters['p3'])
    l4 = linregress(intdate, parameters['p4'])
    l5 = linregress(intdate, parameters['p5'])
    l6 = linregress(intdate, parameters['p6'])
    l7 = linregress(intdate, parameters['p7'])
    l8 = linregress(intdate, parameters['p8'])
    
    # Convert trendline equation to data array
    linearfits = {}
    linearfits['p0'] = []
    linearfits['p1'] = []
    linearfits['p2'] = []
    linearfits['p3'] = []
    linearfits['p4'] = []
    linearfits['p5'] = []
    linearfits['p6'] = []
    linearfits['p7'] = []
    linearfits['p8'] = []
    for thing in intdate:
        linearfits['p0'].append(Linearfnc(thing, l0[0], l0[1]))
        linearfits['p1'].append(Linearfnc(thing, l1[0], l1[1]))
        linearfits['p2'].append(Linearfnc(thing, l2[0], l2[1]))
        linearfits['p3'].append(Linearfnc(thing, l3[0], l3[1]))
        linearfits['p4'].append(Linearfnc(thing, l4[0], l4[1]))
        linearfits['p5'].append(Linearfnc(thing, l5[0], l5[1]))
        linearfits['p6'].append(Linearfnc(thing, l6[0], l6[1]))
        linearfits['p7'].append(Linearfnc(thing, l7[0], l7[1]))
        linearfits['p8'].append(Linearfnc(thing, l8[0], l8[1]))
    
    ## ____ Plots ____ ##
    box = [0.14, 0.14, 0.76, 0.76]
    
    markers = ['o', '^', 'D', 'o', 'D', 's', 'p', 'h', 'd'] 
    colors = ['b', 'g', 'r', 'r', 'g', 'b', 'r', 'g', 'b']
    ylabels = ['A6 Peak Amplitude (in.oz)', 'A6 Peak Location (cm)', 'Gaussian Width (cm)', 
    'A7 Peak Amplitude (in.oz)', 'Gaussian Width (cm)', 'Polynomial Offset (cm)', 'Leading Coefficient', 
    'First Power Coefficient', 'Constant Coefficient']
    titles = ['A6 Peak Amplitude Versus Date', 'A6 Peak Location Versus Date', 
    'Gaussian Width Versus Date for Peak at Bend A6', 'A7 Peak Amplitude (in.oz)',
    'Gaussian Width Versus Date for Peak at Bend A7', 'Polynomial Position Offset Versus Date', 
    'Second Order Polynomial Leading Coefficient Versus Date',
    "Second Order Polynomial's First Power Coefficient Versus Date", 
    "Second Order Polynomial's Constant Coefficient"]
    
    for p in parameters:
        fig = PLT.figure(figsize = (15, 8), dpi = 150)
        ax = fig.add_axes(box)
        ax.set_ylabel(ylabels[int(list(p)[-1])], fontsize = 18)
        ax.set_xlabel('Date', fontsize = 18)
        PLT.title(titles[int(list(p)[-1])], fontsize = 20, y = 1.02)
        ax.grid()
        
        ax.scatter(date, parameters[p], color = colors[int(list(p)[-1])], marker = markers[int(list(p)[-1])])
        ax.plot(date, linearfits[p], '-', color = colors[int(list(p)[-1])])
        fig.savefig(outdir + p + '_vs_date.png')
    
    fig = PLT.figure(figsize = (15, 8), dpi = 150)
    axa = fig.add_axes(box)
    axa.set_ylabel('Gaussian Width (cm)', fontsize = 18)
    axa.set_xlabel('Date', fontsize = 18)
    PLT.title('Gaussian Width Versus Date', fontsize = 20, y = 1.02)
    axa.grid()
    axa.tick_params(labelsize=18)
    axa.scatter(date, parameters['p2'], color = 'r', marker = 'D', label = 'Peak at Bend A6')
    axa.plot(date, linearfits['p2'], '-', color = 'r')
    axa.scatter(date, parameters['p4'], color = 'g', marker = 'D', label = 'Peak at Bend A7')
    axa.plot(date, linearfits['p4'], '-', color = 'g')
    h, l = axa.get_legend_handles_labels()
    #for i, txt in enumerate(filelist):
        #axa.annotate(txt, (date[i], parameters['p4'][i]), position = (date[i], parameters['p4'][i]), size = 'x-small', alpha = 0.5, rotation = 45)
    PLT.legend(h, l, 'upper left')
    fig.savefig(outdir + 'p2_p4_vs_date.png')
    PLT.close(fig)
    
    figb = PLT.figure(figsize = (15, 8), dpi = 150)
    axb = figb.add_axes(box)
    axb.set_ylabel('Amplitude (in.oz)', fontsize = 18)
    axb.set_xlabel('Date', fontsize = 18)
    PLT.title('Amplitude of Peaks at A6 and A7 Versus Date', fontsize = 20, y = 1.02)
    axb.grid()
    axb.tick_params(labelsize=18)
    axb.scatter(date, parameters['p0'], color = 'b', marker = 'o', label = 'A6 Peak')
    axb.plot(date, linearfits['p0'], '-', color = 'b')
    axb.scatter(date, parameters['p3'], color = 'r', marker = 'o', label = 'A7 Peak')
    axb.plot(date, linearfits['p3'], '-', color = 'r')
    h, l = axb.get_legend_handles_labels()
    PLT.legend(h, l, 'upper left')
    figb.savefig(outdir + 'p0_p3_vs_date.png')
    PLT.close(figb)

# ____ Plot fit ___#
def PolyTwoGauss(p, x):           
        # offset = distance between peaks
        # p[0], p[1], p[2] = amplitude , center_x, sigma: for peak A
        # p[3], p[4] = amplitude, sigma: for peak B
        # p[5] = polynomial offset
        # p[6], p[7], p[8] = ax^2 + bx + c
    offset = 52
    y = numpy.zeros(len(x[0]))
    for i, j in enumerate(x[0]):
        y[i] = p[0]*math.exp(-0.5*((j-p[1])/p[2])**2.) + p[3]*math.exp(-0.5*((j-(p[1]+offset))/p[4])**2.) + p[6]*(j - p[5])**2 + p[7]*(j - p[5]) + p[8]
        if y[i] <= 1e-5: 
            y[i] = 1e-5
    return y

def PlotFit(run, outdir):
    """Plot torque vs position data and Gaussian fit to data.
    
    Parameters:
        run (CalTorque): single CalTorque object to plot
        outdir (str): path to folder where plot will be dropped
    """
    date = run.Datetime().strftime('%Y%m%d')
    direction = run.Direction()
    source = run.Source()
    new_parms = run.NewParms()
    position = run.NewPosition()
    torque = run.Torque()
    station = run.Station()
    xdata, ydata = run.BinData()
    
    if station == 'S2':
        xbin = xdata[140:170]
        ybin = ydata[140:170]
    else:
        xbin = xdata[140:-10]
        ybin = ydata[140:-10]
    
    # create folder for plots if one does not exist
    if not os.path.exists(outdir):
        os.makedirs(outdir)
    
    fig = PLT.figure(figsize=(15,8),dpi=150)
    box = [0.14, 0.14, 0.76, 0.76]
    ax1 = fig.add_axes(box)
    ax1.set_ylabel('Torque (in.oz)', fontsize = 18)
    ax1.set_xlabel('Position (cm)', fontsize = 18)
    maxy = max(torque)
    maxx = max(position)
    miny = 0
    
    ## Add station locations
    if maxx >= 328.9:
        ax1.axvline(338.9, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S2', xy = (340.9, miny * 0.8), color = 'r')
    if maxx >= 382.6:
        ax1.axvline(392.6, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S5', xy = (394.6, miny * 0.8), color = 'r')
    if maxx >= 435.7:
        ax1.axvline(445.7, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S8', xy = (447.7, miny * 0.8), color = 'r')
    if maxx >= 491.7:
        ax1.axvline(501.7, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S11', xy = (503.7, miny * 0.8), color = 'r')
    if maxx >= 595.8:
        ax1.axvline(605.8, color = 'r', linewidth = 1, linestyle = '-')
        ax1.annotate('S17', xy = (607.8, miny * 0.8), color = 'r')
    
    ## Add bend locations
    if maxx >= 47.8:
        ax1.axvspan(47.8, 93.5, color = 'g', alpha = 0.2)
        ax1.annotate('Bellows', xy = (48, miny), color = 'g')
    if maxx >= 136.4:
        ax1.axvspan(136.4, 144.6, color = 'g', alpha = 0.2)
        ax1.annotate('First outer bend', xy = (136, miny), color = 'g')
    if maxx >= 196.5:
        ax1.axvspan(196.5, 218.0, color = 'g', alpha = 0.2)
        ax1.annotate('The "S"', xy = (212, miny), color = 'g')
    if maxx >= 228.3:
        ax1.axvspan(228.3, 241.5, color = 'g', alpha = 0.2)
    if maxx >= 274.6:
        ax1.axvspan(274.6, 284.6, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A5', xy = (274, miny), color = 'g')
    if maxx >= 309.0:
        ax1.axvspan(309.0, 318.9, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A6', xy = (309, miny), color = 'g')
    if maxx >= 358.8:
        ax1.axvspan(358.8, 370.8, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A7', xy = (359, miny), color = 'g')
    if maxx >= 414.5:
        ax1.axvspan(414.5, 426.5, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A8', xy = (414, miny), color = 'g')
    if maxx >= 445.3:
        ax1.axvspan(445.3, 455.3, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A9', xy = (445, miny), color = 'g')
    if maxx >= 467.5:
        ax1.axvspan(467.5, 479.5, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A10', xy = (467, miny), color = 'g')
    if maxx >= 523.8:
        ax1.axvspan(523.8, 535.8, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A11', xy = (523, miny), color = 'g')
    if maxx >= 571.6:
        ax1.axvspan(571.6, 583.6, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A12', xy = (571, miny), color = 'g')
    if maxx >= 627.9:
        ax1.axvspan(627.9, 639.9, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A13', xy = (627, miny), color = 'g')
    if maxx >= 652.7:
        ax1.axvspan(652.7, 660.7, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A14', xy = (652, miny), color = 'g')
    if maxx >= 678.5:
        ax1.axvspan(678.5, 681.1, color = 'g', alpha = 0.2)
        ax1.annotate('Bend A15', xy = (678, miny), color = 'g')
    
    PLT.title(date + ' Source ' + source + ' ' + direction + ' Offset Fix by Model Method', fontsize = 20, y = 1.02)
    ax1.grid()
    ax1.tick_params(labelsize=18)
    ax1.plot(xdata, ydata, '-', color = 'maroon', linewidth = 1, label = 'Binned Original Data')
    ax1.plot(xbin, PolyTwoGauss(new_parms, [xbin]), '-', color = 'g', linewidth = 1, label = 'Poly-Two-Gauss')
    #ax1.plot(run.NewPosition(), torque, '-', color = 'blue', linewidth = 1, label = 'Corrected Offset')
    h, l = ax1.get_legend_handles_labels()
    PLT.legend(h, l, 'upper left')
    fig.savefig(outdir + date + '_' + direction + '_' + source + '_modelfit' + '.png')
    PLT.close(fig)
    
def Bin(position, torque, dl = 1, xmin = 0, xmax = 400): # needs to be fixed. ends at 398
    """Bin torque and position data 
    
    Parameters:
        position (array, float): position data
        torque (array, float): torque data
        dl (int): width of bin
        xmin (int): minimum or initial position
        xmax (int): maximum or final position
    
    Returns:
        binned position and torque arrays
    
    """
    xbin = numpy.arange(xmin,xmax,dl)
    ybin = numpy.zeros(xbin.size)
    
    # find index of bin that each position belongs to
    indices = numpy.digitize(position, xbin)
    # Note: indices, xbin, and ybin are parallel arrays
    
    # sort torque values into appropriate bin, max torque as bin value
    for i, t in enumerate(torque):
        if indices[i] < len(xbin) and indices[i] > 0 and abs(torque[i]) > abs(ybin[indices[i]-1]):
            ybin[indices[i] - 1] = torque[i]
    
    return xbin[:-1], ybin[:-1]
    
def frictionModel(x, parms):
    """Friction model
    
    Parameters:
        A (float): constant
        B (float): constant
        mu (float): coefficient of friction
        theta (float): total bend angle
    """
    #import bigfloat
    return parms[0] + (parms[1]*numpy.exp(parms[2]*x))
    
def residual(parms, x, y):
    model_y = frictionModel(x,parms)
    return numpy.sum((model_y - y)**2)
    
def PlotAvgTorquevsAngle(runs, srcTube, dl = 5, maxd = 400, outdir = ''):
    """Plot average torque versus total bend angle
    
    Parameters:
        runs (object): CalTorque run data object to plot
        srcTube (dict): SourceTube x, y data
        dl (int): width of bin in degrees
        outdir (str): path to folder where plot will be dropped
    """
    maxdIndex = (numpy.abs(srcTube['x'] - maxd)).argmin()
    maxB = srcTube['y'][maxdIndex]
    
    binBendAngle = numpy.arange(600,1500,dl)
    binTorque = numpy.zeros(binBendAngle.size)
    allTor = []
    for i in binTorque:
        allTor.append([])
    numTor = numpy.zeros(binBendAngle.size)
    
    # get indices of src tube values in bend angle bins
    # Note: indices and srcTube['x'] are parallel arrays
    indices = numpy.digitize(srcTube['y'], binBendAngle)
    
    # compare position data to bend angle position values
    # to sort torque into appropriate bins
    for r in runs:
        if runs[r].Direction() != 'in': continue #Skip all retractions
        if runs[r].Source() != 'weakTh-228': continue #Skip all other sources
        if 'S5' not in runs[r].Station(): continue # Skip calibration runs

        binTemp = numpy.zeros(binBendAngle.size)
        for i, p in enumerate(runs[r].NewPosition()):
            # find index ot closest position
            index = (numpy.abs(srcTube['x'] - p)).argmin()
            # skip out of bound values
            if p < 0 or p > max(srcTube['x']): continue
            # if closest position is smaller, put in previous bin
            if srcTube['x'][index] < p:
                index -= 1
            # put torque in temporary bin array
            if runs[r].Torque()[i] > binTemp[indices[index]]:
                binTemp[indices[index]] = runs[r].Torque()[i]
        # add temporary bin array to average torque bin array
        for i, b in enumerate(binTemp):
            if b != 0:
                allTor[i].append(b)
                binTorque[i] += b
                numTor[i] += 1
    
    # compute average torque and standard deviation
    avgTorque = []
    bendAngle = []
    maxT = []
    minT = []
    stanDev = []
    for i, t in enumerate(numTor):
        if t != 0 and binBendAngle[i] < maxB:
            avgTorque.append(binTorque[i] / t)
            bendAngle.append(binBendAngle[i])
            stanDev.append(numpy.std(allTor[i], dtype = numpy.float64))
    
    fitParms = [12, 0.00000005, 0.018]
    fit_bounds = [(None,None), (None, None), (None, None)]
    result = minimize(residual,fitParms, args = (numpy.array(bendAngle), numpy.array(avgTorque)), method='TNC', bounds=None,tol=0.0000005,options={'disp':False})
    params = result.x
    
    model = []
    for b in bendAngle:
        model.append(frictionModel(b, params))
    ## TODO convert to radians
    # plot
    fig = PLT.figure(figsize=(15,8),dpi=150)
    ax1 = fig.add_axes([0.14, 0.14, 0.76, 0.76])
    ax1.set_ylabel('Torque (in.oz)', fontsize = 18)
    ax1.set_xlabel('Bend Angle', fontsize = 18)
    PLT.title('Average Torque vs Bend Angle for all weak Th-228 Deployments', fontsize = 20, y = 1.02)
    miny = min(avgTorque)
    maxx = max(bendAngle)
    # add bend locations
    if maxx >= 675:
        ax1.axvspan(675, 720, color = 'r', alpha = 0.2)
        ax1.annotate('Bellows', xy = (0, miny - 2), color = 'r')
    if maxx >= 675:
        ax1.axvspan(679, 760, color = 'b', alpha = 0.2)
        ax1.annotate('First outer bend', xy = (679, miny), color = 'b')
    if maxx >= 711:
        ax1.axvspan(711, 925, color = 'r', alpha = 0.2)
        ax1.annotate('The "S"', xy = (711, miny - 2), color = 'r')
    if maxx >= 915:
        ax1.axvspan(915, 1048, color = 'r', alpha = 0.2)
    if maxx >= 1016:
        ax1.axvspan(1016, 1081, color = 'b', alpha = 0.2)
        ax1.annotate('Bend A5', xy = (1016, miny), color = 'b')
    if maxx >= 1058:
        ax1.axvspan(1058, 1123, color = 'r', alpha = 0.2)
        ax1.annotate('Bend A6', xy = (1058, miny - 2), color = 'r')
    if maxx >= 1085:
        ax1.axvspan(1085, 1163, color = 'b', alpha = 0.2)
        ax1.annotate('Bend A7', xy = (1085, miny), color = 'b')
    if maxx >= 1122:
        ax1.axvspan(1122, 1200, color = 'r', alpha = 0.2)
        ax1.annotate('Bend A8', xy = (1122, miny - 2), color = 'r')
    if maxx >= 1182:
        ax1.axvspan(1182, 1280, color = 'b', alpha = 0.2)
        ax1.annotate('Bend A9', xy = (1182, miny), color = 'b')
    if maxx >= 1269:
        ax1.axvspan(1269, 1348, color = 'r', alpha = 0.2)
        ax1.annotate('Bend A10', xy = (1269, miny - 2), color = 'r')
    if maxx >= 1306:
        ax1.axvspan(1306, 1383, color = 'b', alpha = 0.2)
        ax1.annotate('Bend A11', xy = (1306, miny), color = 'b')
    if maxx >= 1349:
        ax1.axvspan(1349, 1427, color = 'r', alpha = 0.2)
        ax1.annotate('Bend A12', xy = (1349, miny - 2), color = 'r')
    PLT.errorbar(bendAngle, avgTorque, xerr = None, yerr = stanDev, linewidth = 2, elinewidth = 1, capsize=0, ecolor = (.4, .4, 1))
    PLT.plot(bendAngle, model, color = 'black', linewidth = 2, label = 'Fit: y = ' + str(params[0]) + ' + ' + str(params[1]) + 'e^(' + str(params[2]) + 'x)')
    ax1.grid()
    ax1.tick_params(labelsize=18)
    h, l = ax1.get_legend_handles_labels()
    PLT.legend(h, l, 'upper left')
    fig.savefig(outdir + 'AverageTorqueVSAngle_sd.png')
    PLT.close(fig)
